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Section I 
Introduction 


This report describes the activities and accomplishments conducted under contract PO-L- 
14635 with NASA Langley Research Center. Subject matter is the development of a high 
accuracy angular measurement capability for hypersonic wind tunnel facilities located at NASA 
Langley Research Center (LaRC). Specifically, utilization of micro-electro-mechanical sensors 
(MEMS) including accelerometers and gyros, coupled with software driven data acquisition 
hardware, integrated within a prototype measurement system, is addressed. Development 
methodology must address basic design requirements formulated from wind tunnel facility 
constraints and current operating procedures, as well as engineering and scientific test 
objectives. 1 ’ 2 An important component of the development methodology should include a 
description of the analytical framework governing any relationships between time dependent 
multi-axis acceleration and angular rate sensor data and the desired three dimensional Eulerian 
angular state of the test model. 3 ’ 4 A calibration procedure for identifying and estimating critical 
parameters in the sensor hardware is also a key component of the development methodology. 5 ' 7 

Modern experimental and test activities at LaRC demand innovative and adaptable 
procedures to maximize data content and quality while working within a severely constrained 
budgetary and facility resource environment. References 8-9 describe some of these modern test 
procedure philosophies. Unfortunately, many of the LaRC experimental facilities utilize dated 
methodologies and technologies during test activities. Current angular measurement capability in 
LaRC high-speed wind tunnel facilities is a prime example. These current systems are based on 
primitive electro-mechanical measurements of the model drive system and mathematical 
descriptions of flexural characteristics of the model sting. Resulting implications are low quality 
information on the test article pitch, yaw and roll states. In many cases, gross uncertainties enter 
the angular measurement process, which in turn, requires repeated testing and data averaging 
techniques to accomplish test objectives. 
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To overcome these deficiencies, this project focuses on developing a modern and adaptable 
sensor package capable of measuring pitch, yaw and roll while also meeting the size and 
accuracy requirements stipulated by LaRC high-speed short duration wind tunnels. In a general 
setting, pitch, yaw and roll orientation can be indirectly measured from two accelerometer 
packages that sense gravitational acceleration components in various reference frames combined 
with knowledge of the gravity vector orientation relative to the wind tunnel reference system. In 
applications with orthogonality existing between the wind tunnel base and gravity, and with yaw 
rotation about the vertical direction, yaw orientation is unobservable. However, if this yaw 
rotation is about a nonvertical axis, then yaw can be fully recovered. Further, angular orientation 
in all axes can be directly measured by a single gyro based sensor package. Mechanization of 
angular measurements in this project are based on a combination of both approaches. 
Commercially available MEMS sensor technology is utilized to maximize the overall sensor 
package applicability and adaptability to a wide spectrum of wind tunnel facilities. For example, 
use of small MEMS will facilitate sensor placement within the test article itself, rather than at the 
base of the sting support system. In a high-speed facility with significant sting loading, this 
arrangement has an advantage of directly measuring the model motion. Further, as the low cost 
MEMS age or fail, they can be easily replaced with upgraded, functioning units. 

This work is heavily dependent upon NASA LaRC resources, and in particular resources 
allocated to the Models and Systems Branch. These resources include an array of MEMS 
accelerometers and gyros, data acquisition microprocessor cards, dedicated personal computer 
platform, commercial software tools, instrumentation, cables and connectors, high fidelity 
angular calibration tables, and laboratory space to conduct research activities. The project work 
is a research and development effort since the focus is to create a new kind of wind tunnel 
angular measurement system for a unique high-speed tunnel application. The work is also a 
mixed experimental-analytical effort. The experimental activities involve sensor utilization, 
electronic microprocessor data acquisition, software development and coding, and system 
integration. The analytical activities include derivation of required mathematical relationships for 
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multi-axis rotational motion, data processing and extraction, parameterized sensor models, and 
numerical algorithms. Access to NASA LaRC facilities proved challenging for the graduate 
student research assistant, who required escort restrictions at all times. Post September 1 1, 2001 
security issues further slowed access to NASA LaRC. 

The contract Statement of Work (see Appendix A) consisted of fifteen sequential, detailed 
tasks contributing to the overall objective. For purposes of final reporting, these tasks are 
translated into the following areas. 

1 . Survey tunnel angular measurement requirements and available MEMS performance. 

2. Derive analytical relationships for angular motion and sensor characteristics. 

3. Develop accelerometer and gyro sensor calibration and attitude construction procedures. 

4. Design data acquisition hardware-software and conduct system integration and test. 

These areas are briefly outlined here before moving on to the dedicated sections with detailed 
reporting of the activities. Emphasis during the first two thirds of the project was given to 
achieving delivery of a working prototype system at contract closure. After discussions with 
NASA LaRC, the final third of the project was focused on gyro calibration procedures. All key 
aspects of the angular measurement system prototype have been investigated. Due to limited 
commercially available or LaRC owned MEMS sensor performance, a fully working system that 
meets all design performance objectives was not achieved. However, the prototype system was 
sufficiently explored such that, if given appropriate MEMS performance, the prototype system 
would in fact satisfy the objectives. At contract closure, the prototype system consists of loosely 
integrated hardware and software components, many of which are NASA LaRC property and are 
thus not part of the deliverables. The main deliverables of the contract are the data acquisition 
software, sensor calibration software, attitude construction software, and final report describing 
how these components are integrated into an overall working system. 

Section II describes the first area dealing with design requirements and sensor performance. 
High-speed wind tunnel operational constraints and environmental factors were collected for 
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purposes of bounding the measurement system design space and characterizing the current 
angular measurement fidelity. Wind tunnel data is based on a single but representative LaRC 
high-speed facility, the 31 -Inch Continuous Flow Hypersonic Mach 10 Tunnel. This data was 
collected from a detailed facility tour coupled with in-depth discussions with technician staff. 
Constraint information covered operational, temporal, thermal, pressure, dimensional, electronic, 
range, and accuracy aspects. Baseline angular measurement accuracy is stated to be ±0. 1 deg in 
pitch and ±0.5 deg in yaw. However, these estimates are most likely optimistic, with more 
realistic values being near ±1.0 deg in pitch and yaw. Desired resolution levels for high-speed 
aerodynamic flow phenomenon were also collected and estimated to be near ±0. 1 deg. Based on 
vendor specification sheets, a survey of commercially available MEMS accelerometer and gyro 
units was also conducted. Some of these units are within the NASA LaRC hardware inventory. 
Based on manufacturer’s stated performance, a first order estimate indicates the prototype 
MEMS based measurement system using accelerometers should attain an approximate accuracy 
of ±0.1 deg in pitch, yaw and roll. Similarly, a prototype system using gyros (assuming stable 
sensor characteristics, high fidelity calibration, and efficient pretest activities) is predicted to 
achieve approximately ±0.1 deg accuracy in pitch, yaw and roll. Unfortunately, as discussed 
latter, performance of the available MEMS gyro units fell significantly short of the 
manufacturer’s claims. 

Section III describes the second area dealing with the underlying analytical framework. In 
this section, appropriate reference frames for important components of the wind tunnel and 
angular measurement system are introduced. An Eulerian description of the test article attitude 
state, consistent with the tunnel hardware, is defined. Because of the limited hardware range of 
rotation, mathematical singularity conditions lie outside the tunnel operational limits. Thus, an 
Eulerian description for the angular state is appropriate in this application. A general rotation 
sequence, which is consistent with the 31 -Inch Hypersonic Tunnel, is selected to facilitate multi- 
axis pitch, yaw and roll attitude estimation. Accelerometer based measurements will be used in 
an indirect and sequential procedure for extracting attitude information by sensing the 
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gravitational acceleration components first in the sting axes and second in the model axes and 
comparing with the known gravity vector orientation. Therefore, governing relationships 
between accelerometer sensed gravity components and the model angular state are developed 
and documented. In contrast, gyro based measurements will provide a direct procedure for 
extracting attitude information by sensing attitude rates and conducting integration. Fundamental 
equations describing gyro sensed multi-axis angular rates and their relationship to the model 
angular state are generated. Accelerometer and gyro sensor characteristics are also 
mathematically parameterized in this section. Accelerometer and gyro parameter models 
considered in the project include mechanical-electrical hardware aspects such as bias, sensitivity, 
and azimuth and coning misalignment angles. These parameters are considered in all axes. The 
analytical relationships derived in this section provide the basis for pre-measurement sensor 
calibration activities, as well as post-measurement attitude construction efforts. 

Section IV describes the third area dealing with development of sensor calibration and 
attitude construction procedures and techniques. Development of calibration and construction 
algorithms was a major effort within the overall project. Gravity sensing accelerometer 
calibration, as well as angular rate sensing gyro calibration, can be mathematically described by 
a nonlinear algebraic relationship. Recommended sensor calibration procedures are based on an 
existing technique thoroughly described in Reference 5. This technique uses an iterative Gauss- 
Newton numeric solution strategy and was implemented by the contractors using Math Works 
MATLAB. 10 These calibration procedures were tested with simulated data and found to be 
reliable and accurate. Gravity sensing accelerometer based attitude construction is also described 
mathematically with a nonlinear algebraic relationship. Here, an analytic closed-form solution 
for the attitude construction is possible, as well as a Gauss-Newton iterative solution. Both 
solutions were implemented by the contractors using MATLAB. Rate gyro based attitude 
construction is governed by a set of nonlinear differential relationships, thus making the attitude 
construction process here fundamentally distinct from the accelerometer based construction task. 
A direct Runge-Kutta numerical integration strategy was found to be effective in solving the 
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attitude construction equations. This technique was also implemented by the contractors in 
MATLAB. These attitude construction procedures were again tested with simulated data and 
found to be reliable and accurate. Note these full-blown calibration and construction procedures 
were not tested in the laboratory with real data and within the data acquisition control software 
environment. However, simplified single-axis versions of these concepts were tested in the 
laboratory with real data. 

Section V documents the final area dealing with design of data acquisition hardware- 
software and system integration and test. For the prototype system, no freedom existed in 
hardware design or selection. Utilization of NASA LaRC supplied hardware was mandated. 
Major hardware components consisted of a Pentium II class personal computer with video 
monitor, National Instruments 4 channel 200 kHz NI-4452 data acquisition card, Analog Devices 
ADXL105 MEMS micromachined polysilicon vibrating structure accelerometer, and muRata 
ENC-03J piezoelectric vibrating gyroscope or rate gyro. A wide range of design freedom and 
selection existed with data acquisition controlling software. National Instruments LabVIEWU.12 
was selected for this function based on its powerful capabilities, ease of use, compatibility with 
the NI-4452 data card, and acceptance at NASA LaRC. LabVIEW’s graphical user interface was 
used to create multi-functional control instrumentation that reads input signals from the data 
card, processes and filters the input signals, calibrates several key sensor parameters, transforms 
the acceleration and angular rate data to angular position states, records data to storage devices, 
and outputs real-time information to a video display in numeric and graphical formats. Real-time 
output can be monitored by the wind tunnel technicians and experimental aerodynamicists to 
evaluate ongoing test conditions. Software design and testing was a major effort within the 
overall project. A significant problem that was encountered related to critical LabVIEW 
functions interfacing with the NI-4452 card. To circumvent this problem, several core functions 
from Wyle Laboratories were utilized in the control software. Integration of the hardware and 
software components to a working end-to-end simplified angular measurement system prototype 
was successfully accomplished. This end-to-end operational capability was tested and 
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demonstrated in the laboratory. Note only simplified single-axis calibration-construction 
capability was demonstrated here. Angular measurement accuracy objectives using available 
MEMS accelerometers was achieved. Unfortunately, these objectives were not achieved with 
available MEMS gyros. Large performance discrepancies between manufacturer published data 
and observed data were discovered. A major conclusion of the work is that better performing 
gyros are required. 
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Section II 

Design Requirements and Sensor Performance 

A. Wind Tunnel Constraints 

Design of a prototype wind tunnel Angular Measurement System (AMS) depends strongly 
upon the characteristics and constraints imposed by the wind tunnel facility itself. Further, such 
information will benchmark the current AMS accuracy for the purpose of latter comparison with 
the new AMS system. In this project, the 31 -Inch Continuous Flow Hypersonic Mach 10 Tunnel 
is selected as a basis for the AMS design because this facility is frequently used and is 
representative of several other LaRC high-speed wind tunnel facilities. Wind tunnel constraints 
and characteristics were collected from a detailed facility tour coupled with in-depth discussions 
with technician staff, as well as internet information on the tunnel. Detailed information was 
collected and recorded in several areas including operational, accuracy, dimensional, temporal, 
dynamic, range, thermal, and pressure aspects. 

Figure 1 illustrates the sequence of events during a test in the 31-Inch Hypersonic Tunnel. 
The aerodynamic model is prepared for testing while the injection housing is in the exposed 
position. The test article is mounted on the pitch-yaw sting drive system and the system is rotated 
to the desired initial angles. Fixed roll orientation is achieved by direct sting mounting. Figure 2 
shows the sting drive system arrangement. Instrumentation attached to the model is verified for 
correct and stable readings. Electro-mechanical measurements are taken of the model orientation. 
The injection housing is then rotated to the closed position while the wind tunnel flow conditions 
are adjusted to the desired values. When the flow conditions are met, the model is injected 
rapidly into the high-speed flow. As the model is injected into the test chamber, high-speed 
airflow loads the drive system and flexes the sting. Consequently, the model orientation during 
test is different from the measured pre-injection condition. Correction factors based on models of 
the drive system hysteresis and sting flexural characteristics are commonly used in an attempt to 
improve knowledge of the model attitude test conditions. However, significant uncertainties 
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remain in the baseline measurement process. Baseline angular measurement accuracy was stated 
to be ±0.1 deg in pitch and ±0.5 deg in yaw. Further discussion revealed these estimates are most 
likely optimistic, with more realistic values being near ±1.0 deg in pitch and yaw. More 
important than the exact value for angular accuracy is the knowledge that sizeable unknowns 
exist in the current AMS. The new AMS should provide an accuracy level significantly better 
than existing system performance described here. Thus, a requirement for reliable angular 
measurement accuracy of ±0.1 deg in all axes will be considered for the AMS design. 










Now consider dimensional constraints placed on the AMS design. Typical model size in 
the 31 -Inch Hypersonic Tunnel is 12-15 in. Models are constructed to meet center of gravity and 
moment of inertia requirements precisely. Sting mounting is accomplished with rear attachment 
geometry as shown in Figure 2 using several standardized interface dimensions. To circumvent 
any necessary altering of the model and its associated inertial characteristics, and to avoid 
variability with each model, mounting of the sensor package internally within the sting tip is 
recommended. Such mounting will minimize modification to the existing facility and test 
procedures, and in particular interfacing with the test article. Stings for the 31 -Inch Hypersonic 
Tunnel are typically slender, hollow cylindrical 1-1.5 ft long tubes constructed from various 
metallic alloys. Internal sting diameters range between 0.5 and 0.75 in. Volume and dimension 
constraints placed on the sensor package are thus quite small, and the weight is also suggested to 
be as light as possible. An ability to reliably secure and unsecure the package within the sting tip 
is also desired. Electrical wiring to and from the sensors can also be accomplished through the 
hollow sting shaft. 

Next, consider constraints from test duration and drive system angular range. Typical test 
times for the 31 -Inch Hypersonic Tunnel will be less than 2 min. This duration does not pose any 
significant constraint on the AMS design. However, the model preparation time during which the 
injection housing is in the exposed position can be significantly longer than the actual test time. 
Preparation time could be an order of magnitude longer. Since final alignment of the AMS will 
be performed somewhere in this window, further constraints are placed on the AMS design. The 
drift characteristics of the MEMS accelerometer and gyro units, in combination with drift 
estimation accuracy, must accommodate durations on the order of 10-15 min in meeting the 
angular measurement requirements. The pitch-yaw drive system has an operational rotation 
range of ±60 deg in pitch and ±15 deg in yaw. Common testing angles, however, fall within ±15 
deg pitch and ±5 deg yaw. Roll capability is unlimited (±180 deg), but again is commonly held 
within ±15 deg. It is important to note that pitch and yaw orientation can be varied during the test 
while roll is fixed. The model is usually held for 2-3 s at each test point with transition time 
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between test points requiring approximately the same duration. This dynamic motion places a 
mild constraint on the required dynamic response of the MEMS units and data acquisition 
sample time. Modern sensors and data input cards will have no problem in satisfying these 
constraints. Because of the small hardware range of rotation, there are no constraints for using a 
singularity free, non-minimal angular state description. Thus, a simpler Eulerian description for 
the angular state will be acceptable. 

Finally, the harsh environment of the hypersonic test chamber requires the sensor package 
to be functional in a high temperature, vacuum environment. Temperature within the 31-Inch 
Hypersonic Tunnel test chamber commonly rises to 1,350 °F. This temperature is well beyond 
the operational range of any MEMS gyroscope or accelerometer. However, the expected 
temperature within the sting tip, which itself is inside the test article is significantly lower than 
test chamber temperatures. The actual temperature within the sting tip is not provided but is 
estimated to be 100-300 °F. This temperature is on the fringe of most commercially available 
MEMS temperature range capabilities. Further, MEMS units are typically temperature sensitive. 
Thus, another constraint imposed from the tunnel facility is for thermal compensation of some 
type, and detailed thermal calibration. The test chamber also experiences a vacuum condition on 
the order of 1-2 psi during operation. Constraints for acceptable operational performance in a 
vacuum will also be placed on any AMS design. Due to limited project resources, thermal and 
pressure characteristics in the design of the prototype AMS are not considered, and will have to 
be addressed in follow on efforts. 

B. Aerodynamic Test Objectives 

Design of the prototype wind tunnel AMS also depends strongly on the engineering and 
scientific test objectives. NASA LaRC high-speed wind tunnel facilities are typically used to 
explore and quantify the aerodynamic characteristics of advanced hypersonic concepts, air-to-air 
missiles, orbital launch systems, and re-entry vehicles over a wide range of speeds and 
aerodynamic attitudes. For example, hypersonic wave rider concepts in the cruise configuration 
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would be tested at very low (near zero) angles of attack and side slip angles, while maneuvering 
re-entry vehicles would be tested at large (possibly approaching 45 to 60 deg) angles of attack 
and side slip angles. Associated flow physics are complex and may exhibit attached oblique 
shocks, detached bow shocks, large regions of separated chaotic flow, body generated vortices, 
and thick turbulent boundary layers. These characteristics may be highly sensitive to and 
dependent on the aerodynamic angles, thus requiring a high fidelity model orientation 
measurement capability. Reference 2 implies current industry-government experimental program 
goals are to achieve an angular wind tunnel measurement resolution of ±0.1 deg or better. 
Therefore, an engineering and scientific test imposed requirement for reliable angular 
measurement accuracy of ±0.1 deg in all axes will be considered for the AMS design. This 
requirement is consistent with previous accuracy requirements driven by current wind tunnel 
facility capabilities and the need to improve upon this capabilities. 

C. MEMS Survey 

Several MEMS accelerometers under consideration for the prototype AMS are listed in 
Table 1. Many other MEMS accelerometers are commercially available, but this report will 
focus only on those units listed in Table 1. These units include Analog Devices’ ADXL105, and 
MEMSIC’s MX201A, MXA2500A and MXD2020A units. NASA LaRC has a large supply of 
the ADXL105 units but currently does not have any units from MEMSIC. Detail sensor 
specifications can be found in References 13-14. The ADXL105 MEMS accelerometer uses an 
open-loop architecture consisting of a micromachined polysilicon vibrating structure outfitted 
with a central capacitance plate and positioned between two other fixed capacitor plates. The 
fixed plates are used to nominally vibrate the central plate, as well as to sense unsymmetric 
motion resulting from accelerations. The MX201A, MXA2500A and MXD2020A MEMS 
accelerometers measure internal heat transfer changes within a gaseous proof mass through 
natural convection caused by acceleration. A gaseous cavity is thermally loaded by a heater bar 
and instrumented with several thermocouples. Acceleration results in an unsymmetrical thermal 
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Table 1 

MEMS Accelerometer Specifications 

Vendor 

Analog Devices 

MEMSIC 

MEMSIC 

MEMSIC 

Model 

ADXL105 

MX201A 

MXA2500A 

MXD2020A 

Linear Resolution (mg) 

2 

1 

<2 

<2 

Angular Resolution (deg) 

0. 1 146 

0.0573 

<0.1146 

<0.1146 

Sensitivity (mv/g) 

250 

500 

500 

- 

Max Acceleration (g) 

±5 

±1 

±1 

±1 

Bandwidth (Hz) 

10k 

40 

30 

30 

Noise (pg/Hz 1/2 ) 

225 

750 

750 

40 

Dimension (h/w/d mm) 

10.64/12.07/5.08 

5/5/2 

5/5/2 

5/5/2 

Mass (g) 

<2 

< 1 

< 1 

< 1 

No. of Axes 

1 

2 

2 

2 

Temperature Range (°C) 

0/+70 

-40/+125 

-40/+105 

-40/+105 

Temperature Sensor 

Included 

Included 

Included 

Included 


gradient in the cavity leading to a thermocouple voltage signal that is proportional to 
acceleration. These sensors are also open-loop devices. 

All four accelerometer models are sufficiently small to satisfy the dimensional constraints 
from the internal sting diameter. Even though thermal effects are not considered in the AMS 
design at this time, the MEMSIC units are noted to provide a somewhat larger but still marginal 
window for temperature range. No information is available on acceptable pressure ranges. Linear 
resolution of all four accelerometers range between 1 to 2 mg, and sensitivity ranges between 
250 and 500 mv/g. For the gravity sensing wind tunnel application, this resolution can be 
converted to an angular resolution. Linear resolution is converted to angular resolution through 
an inverse sine operation. In the case of the ADXL105 accelerometer, 2 mg linear resolution 
corresponds to an angular resolution of 0. 1 146 deg (A0 = sin~*(R + sin(©)) - © with R = 0.002 
and 0=0 deg). Note the angular resolution is a function of geometry and is best when the 
sensing axis of the accelerometer is orthogonal to the gravity vector. When the accelerometer is 
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rotated, the angular resolution degrades as the transformation from linear to angular occurs on a 
flatter region of the sine function curve. For ADXL105, the resolution becomes 0.1 187 deg when 
the accelerometer is rotated by 15 deg. With multi-axis sensing, resolution of any of the four 
candidate accelerometer units is very close to meeting the angular measurement requirement of 
±0.1 deg. Dynamic response performance noted in Table 1 should easily satisfy any constraints 
from the drive system duty cycle motion. Also, test duration related constraints are not critical 
here since the accelerometer signals do not require integration for angle construction. 

Next, a set of MEMS gyroscopes listed in Table 2 are considered. Many other MEMS 
gyros are commercially available, but this report will focus only on those units listed in Table 2. 
These units include muRata’s ENC-03J and ENV-05F-03, and Samsung’s MGL80301ZT0 and 
MG-L106C units. NASA LaRC has a large supply of the ENC-03J units but currently does not 
have any units from Samsung. Detail sensor specifications can be found in References 15-16. 
The gyroscopes under consideration are rate gyros which give output signals proportional to 
sensor angular rate. The ENC-03J and ENV-05F-03 MEMS gyros use an open-loop architecture 
consisting of a ceramic bimorph vibrating structure excited electrically. When this oscillating 
unit experiences an angular rate, an electrical signal is generated proportional to the resulting 
Coriolis force. The MGL80301ZT0 and MG-L106C MEMS gyros function in a similar way but 
incorporate proprietary material and manufacturing techniques. These sensors are also open-loop 
devices. 

Table 2 shows the ENC-03J, MGL80301ZT0 and MG-L106C gyros satisfy the 
dimensional constraints. Note the ENV-05F-03 gyro width dimension exceeds the sting internal 
diameter constraint (0.772 in width vs. 0.75 in diameter). Even though thermal effects are not 
considered in the AMS design at this time, both the muRata and Samsung units provide only a 
marginal window for operation within an elevated thermal environment. No information is 
available on acceptable pressure ranges. Although not directly supplied by the manufacturer, an 
angular rate resolution value can be estimated by dividing sensitivity into noise. Only the 
muRata ENV-05F-03 gyro gives a reasonable resolution when compared with expected model 
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Table 2 MEM 

S Gyro Specifications 

Vendor 

muRata 

muRata 

Samsung 

Samsung 

Model 


ENV-05F-03 

MGL80301ZT0 

MG-L106C 

DC Drift (deg/s) 

<5 

<9 

0.2 

0.2 

Sensitivity (mv/deg/s) 


25 

0.6 

0.6 

Max Angular Rate (deg/s) 



±500 

±500 

Bandwidth (Hz) 

50 

7 

15 

50 

Noise (mv) 

- 

5 

18 

18 

Dimension (h/w/d mm) 



8 / 15 / 6.55 

8 / 15 / 6.55 

Mass (g) 

1 

20 

- 

- 

No. of Axes 

1 

1 

1 

1 

Temperature Range (°C) 

-5/+75 

-30/+80 

-5/+55 

• 51+15 


drive rate schedules. The gyros under consideration show significant output signal drift. The 
Samsung gyros exhibit better drift performance than muRata units, but the maximum drift value 
indicates accumulation of 0.2 deg/s drift over a 2 min test results in 24 deg offset. Recall model 
preparation time could lead to significantly larger test durations which only exacerbate the 
problem. If 1) gyro drift characteristics of 0.2 deg/s were highly repeatable and stable, 2) drift 
calibration performance is on the order of 1% (±0.002 deg/s knowledge of actual drift), and 3) 
operational duration could be limited to 120 s with efficient test procedures, then and only then 
does the angular measurement accuracy being to approach the ±0.1 deg requirement. Under these 
drift assumptions, only the ENV-05F-03 gyro is theoretically feasible for satisfying angular 
resolution objectives. Unfortunately, as discussed in Section V-C, the actual drift value observed 
in test for the available gyros was much higher than stated by the manufacturer, and more 
importantly, the drift characteristics were highly variable. Therefore, drift performance for any 
gyro based prototype AMS will be a critical parameter. Finally, dynamic response performance 
noted in Table 2 should easily satisfy any constraints from the drive system duty cycle motion. 
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D. Prototype Angular Measurement System 

A rough schematic of the prototype AMS is described in Figure 3. The architecture in 
Figure 3 shows the hardware and software components, their relationship to one another, and the 
major functions occurring within the AMS prototype. The AMS operates in one of two modes: 
sensor calibration or attitude construction. AMS hardware includes an accelerometer-gyro sensor 
package, data acquisition board, personal computer with video monitor, and electrical 
connections between these devices. Major functions performed by these components include 
motion detection and motion measurement, data sampling and analog-to-digital signal 
conversion, engine for real-time and post-measurement processing, and data storage and data 
display. AMS software resides entirely within the computer and incorporates multiple languages 
with command line and graphical interface modes. Although the prototype AMS can not function 
properly without the presence of high performance hardware components such as motion sensors 
and acquisition devices, as can be observed in Figure 3, software is the most critical component 
in AMS operation. In the sensor calibration and attitude construction modes, software functions 
include data processing and filtering, data transfer, storage and display, estimation of sensor 
characteristics (calibration), and Euler angle reconstruction from raw data. Note the attitude 
construction mode requires information from the sensor calibration mode. Further note the 
attitude construction mode can operate in a lower functionality real-time submode or a higher 
functionality post-test submode. 

Utilization of the AMS is envisioned as follows. Sensor calibration activities will be 
conducted in the Instrumentation Laboratory of the Models and Systems Branch at NASA LaRC. 
This activity will make use of the Dual Axis or Single Axis High Precision Angle Indexing 
Tables in the lab. Once the calibration parameters are identified, the AMS will be transported to 
the appropriate hypersonic wind tunnel facility. The AMS will then be integrated with the test 
article and tunnel data acquisition system to perform attitude construction functions during and 
after the actual test. Independently obtained angles from the accelerometer sensors and gyro 
sensors can be used in a redundancy role to monitor for corrupted data in either channel, or they 
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- Sensor Calibration 



Computer 


Figure 3 Angular Measurement System 


can be fused in such a fashion to supplement either source’s weaknesses, although such functions 
were not considered in this project. At this preliminary development stage, the prototype system 
consists of loosely integrated hardware and software components, many of which are NASA 
LaRC property that are currently being utilized in other projects. At one time during the later 
stages of the contract, these components were in place and formed a working prototype 
demonstration. Only simplified single-axis laboratory calibration and construction capability was 
demonstrated. Currently, the AMS prototype components are disassembled and therefore not 
able to perform working system demonstrations. However, this report describes the working 
system details and performance, and how it could be easily re-assembled. AMS prototype 
characteristics of this sort are consistent with the main project objective of exploring and 
demonstrating the feasibility of an AMS capability, as opposed to constructing a professionally 
packaged and fully completed system. 
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Section III 

Analytical Framework 


A. Angular State Description 

An Eulerian description for the model angular state will be utilized by the AMS. This 
formulation provides advantages by working with insightful physical variables and a minimum 
state description. Singularity conditions are automatically avoided by the limited pitch-yaw drive 
system operational range. Right-handed coordinate systems used in this report are described in 
Figure 4. The tunnel fixed coordinate system (x t y t z t ) is regarded as an inertial frame, while the 
model (body) fixed coordinate system (x m y m z m ) and sting fixed coordinate system (x st y st z st ) are 
rotated relative to the tunnel fixed coordinate system. These systems are oriented conventionally: 
x t along tunnel axis (flow direction), y t orthogonal to tunnel axis (horizontal), z t orthogonal to 
tunnel axis (vertical), x m along model longitudinal axis (nose), y m along model lateral axis (right 
wing), z m along model normal axis (down), x st along sting axis (forward), y st along pitch drive 
axis (horizontal), and z st along yaw drive axis (tilted vertical). For simplification purposes, the 
wind tunnel and associated reference frame are assumed to be oriented relative to the earth fixed 
coordinate system (also assumed inertial) such that the gravitational acceleration vector points 
along the -z t direction. If a more general geometric condition is required, the mathematical 
framework could be easily modified. 

Initial orientation of the sting fixed and model fixed axes is such that x st = x m = -x t , y st = 
y m = y t , and z st = z m = -z t . Note at the initial orientation, axes are aligned but x st = x m and z st = 
z m have opposite directions from x t and z t . This initial sting fixed and model fixed orientation is 
obtained from the tunnel fixed axes by a +180 deg rotation about y t . A non-standard three angle 
Euler rotation sequence is used to describe the transformation relationship between the initial and 
final orientations of the sting fixed and model fixed reference frames. To transform from initial 
orientation to the sting and model axes, the Euler rotation sequence is specified as pitch-yaw-roll 
(0-vp-O). First, the initial orientation is rotated an angle 0 about y t . Second, the resulting 
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intermediate axes are rotated an angle 'P about the corresponding z axis. This rotation sequence 
leads to the sting fixed axes. Third, sting fixed axes are rotated an angle O about the x st axis. The 
indicated Euler rotation sequence leads to the model fixed axes. This rotation sequence is 
consistent with the 31-Inch Hypersonic Tunnel model sting drive system and test procedure. In 
this facility the azimuth (yaw) drive system rides on the elevation (pitch) drive system. Further, 
the test model is attached to the sting with a fixed rotation (roll), which can be interpreted as a 
roll drive system riding on the yaw drive system in an advanced state of the art tunnel facility. 
Roll angle can always be set to a known fixed value to represent specific characteristics of the 
31 -Inch Hypersonic Tunnel. Figure 4 describes the pitch-yaw -roll angle sequence. 



Figure 4 Tunnel-Sting-Model Axes and Euler Sequence Description 


Misalignment between the sting and/or model fixed axes and the sensor fixed axes can 
exist due to attachment and mounting inaccuracies associated with sting-model-sensor 
interfacing, as well as internal sensor manufacturing inaccuracies. Here, sensor fixed axes refer 
to the true internal excitation axes within the sensor. Such misalignment will be described by the 
coning angle (C) and azimuth angle (A). Figure 5 illustrates the definition of these two 
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misalignment angles for the case where x st (sting x axis) is the desired measurement axis and x sn 
(sensor x axis) is the true measurement axis. Coning angle is defined as the angle subtended by 
the desired and true measurement axes. A cone shape is created when the excitation axis is 
rotated about the desired axis. To define the azimuth angle, consider projecting the cone and 
excitation axis onto the y st -z st plane. The projected shape is a circle with a radial marking. 
Azimuth angle is defined as the angle subtended by the -z st axis and projected excitation axis. 
The internal sensor fixed excitation axes relative to the sting fixed axes can be defined by a 
rotational transformation associated with the coning and azimuth angles. This transformation is 
geometrically defined by 1) heading (A) around the circle relative to y st , and 2) spin (C) about 
the heading radial. Other important parameters associated with the sensor behavior are bias (B) 
and sensitivity (S). A systematic deviation exhibited by a measurement value from a reference 
value is defined as bias. Further, sensitivity is defined as the ratio of observed output signal 
produced with respect to a specified unit of sensor input change. 



Figure 5 Azimuth and Coning Angle Description 
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B. Accelerometer Equations 

Derivation of the underlying analytical framework between multi-axis gravitational 
acceleration data measured by two accelerometer sensor packages and the desired three 
dimensional Eulerian angular state of the test model is considered in this subsection. This 
analytical framework provides the basis for accelerometer sensor calibration and attitude 
construction logic presented in Section IV. With the geometry described in Figure 4, the 
gravitational acceleration vector (algebraic form) expressed in the tunnel fixed reference frame is 

G t = [ 0 0 -G] T (1I1-B.1) 

where G denotes the local value of gravity. G t is considered to be the known reference vector. 

General orientation of a three axis set with respect to an orthogonal reference frame is 
defined by a minimum of three angles. However, the general orientation of a vector with respect 
to this same reference is defined with a minimum of two angles. Therefore, extraction of the 
complete pitch-yaw-roll orientation state from a gravity vector measurement scheme employing 
a single accelerometer sensor package is not feasible. For example, as the generally oriented 
model fixed frame is rotated about the gravity vector, an infinite set of pitch-yaw-roll solutions is 
generated corresponding to a single set of model fixed gravity components. Therefore, in a 
general framework, gravitational acceleration will be measured in both the sting fixed and model 
fixed axes using two accelerometer sensor packages. By relating these measured acceleration 
components with known components in Equation (III-B.l) in a two-step process, the full angular 
state can be determined. This approach is applied to development of the general angular 
measurement system. A simplified case, with known roll angle and only pitch-yaw construction 
based on a single sensor package mounted in the sting fixed axes (i.e., the 31 -Inch Hypersonic 
Tunnel application), can always be extracted from the general framework. 

Consider the following sequential transformations of the gravity vector in tunnel fixed axes 
through the model fixed axes. In Equations (III-B.2)-(III-B.5), Gj denotes the gravity vector 
expressed in frame i, G x .-G y .-G z . denote components of Gj, Tj/j denotes the transformation 
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matrix from frame j to frame i, and i, j = t, i, 1 1 , st, m representing tunnel fixed, initial 
orientation, first intermediate, sting fixed, and model fixed axes, respectively. 


G i = T„,G, G i =[G x .G y .G Z| ] T 

G il =T il/i G i G U = [ G \i] G Vj] G z n ] 

G st = T st/il G il G st = G 'st G yst Gz st 

G m = T m/st G st G m = [ G x m G y m Gz m] 


’-10 o' 

T i/t = 0 1 0 (I11-B.2) 

[o 0 -1 

cos(0) 0 -sin(0) 

T il/i = 0 1 0 (1II-B.3) 

sin(0) 0 cos(0) 


cos(^) sin('P) 0 

T st 41 = -sin(W) cos(T') 0 (11I-B.4) 

0 0 1 


1 0 0 

T m/st = 0 cos(O) sin(O) (H1-B.5) 
0 -sin(O) cos(O) 


The gravitational acceleration vector expressed in the sting fixed reference frame in terms of the 
known gravity vector is obtained by combining Equations (III-B.2)-(HI-B.4), or 


G st “ Zst/il^ii/iTjA G t 


(11I-B.6) 


The overall transformation matrix T st / t from tunnel fixed to sting fixed axes in Equation (I1I-B.6) 
is obtained from a product of the component transformation matrices. Note the gravitational 
acceleration vector expressed in the model fixed reference frame as a function of the gravity 
vector expressed in the sting fixed axes is given by Equation (III-B.5). 

Now consider effects from misalignment errors due to the azimuth and coning angles in the 
sting mounted accelerometer sensor package. The following analysis must be conducted 
individually for each desired measurement axis. First consider the sting x st axis as the desired 
measurement axis (Figure 5). A sequential transformation consisting first of azimuth rotation 
(A Xst ) about x st leading to a second intermediate frame, and next of coning rotation (C Xs( ) about 
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the intermediate y axis which leads to the true measurement or sensor fixed reference frame 

(x sn y S n z sn ) associated with x st , is implemented. For the y st axis, azimuth rotation (A y ) is 

about y st and coning rotation (Cy s{ ) is about the intermediate z axis leading to the y st specific 

sensor fixed axes (x sn y sn z sn ). Finally, for the z st axis, azimuth rotation (A z ) is about z st 
yst yst yst st 

and coning rotation (C Zst ) is about the intermediate x axis leading to the z st specific sensor fixed 

axes (x sn y sn z sn ). These transformation relationships are summarized in Equations (III- 

z st z st z st 

B.7)-(III-B.12) where Gj, G x .-G y .-G z ., and Tj/j have the same meaning as before, but with i, j = 
st, i2k, sn k representing sting fixed, second intermediate, and sensor fixed axes, respectively, and 

k = x st- y s t> and z sc 


G 


' x st 


G i2 x " T i2 x ,/st G st 
x st x st 


_ - 

T 

1 0 0 

G G G 

T i2 x /st = 
x st 

0 cos(A % ) sin(A Xst ) 

st st st 

0 -sin(A Xst ) cos(A Xst ) 


(II1-B.7) 


° S "% Ts "x s ,' l2 x st G|2 % t 


(III-B.8) 


G 


s " x st 





l sn x Uly 
x st ^st 


c°s(C Xst ) 0 -sin(C Xst ) 

0 1 0 

sin(C Xst ) 0 cos(C Xst ) 


G 


i2 


yst 



/fet G st 


(III-B.9) 


G 


i2 


y s t 





l i2 v /st “ 

yst 


cos(A Jst )0-sin(A yst ) 

0 1 0 
sin(A >si ) 0 cos(A )st ) 


rsn y s t Tsn ys/ ,2 yst Gi2 yst 


(III-B.10) 
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sn 


yst 





l sn v /i2v 
*st ^st 


cos(C yst ) sin(C yst ) 0 

“ sin(C y s t ) cos(C y s t ) 0 

o 0 1 


Gn 


-z St 


“ T i2 z /st G st 
z st 


(III-B.l 1) 
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0 0 
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‘st 


The gravitational acceleration vector expressed in the three sensor fixed reference frames in 
terms of the sting fixed gravity vector is obtained by combining Equations (III-B.7)-(III-B.12), or 


1 _ rp rp p 

r sn Y x sn Y M x i2 Y /st VJ st 
x st x st x st x st 


r siiy^ T sn yst yi2 ys Ti2 y ^/st^st 


l sn Y /st 
x st 


sn y s t /st 


'sn 


*st 


= T, 


T: 


G c 


sn, 42, A i2, /st^st 


l sn, „/st 
z st 


(III-B.13) 


The overall transformation matrices T sn /st , T sn /st , and T sn /st from sting fixed to sensor 

x st y s t z st 

fixed axes in Equation (III-B.13) are obtained from a product of the component transformation 
matrices. Note only the x, y, and z components, respectively, in these three vector expressions 
are of interest. Extracting these three components and constructing a new vector relationship 
yields Equation (III-B.14). 


^sn st ^sn st /st^st 


Gsn st 


G x G G 

sn st >sn st sn st 


sn st /st 


sn x st *ti 

r S n y st /fet 2 

r9l z st /st 3 


(III-B.14) 


In Equation (III-B.14), T s „ x T sn ^ / st2 , and T SI , z /st3 denote the first, second, and third rows 
of T sn Xst /st’ T sn Vst /st’ and T sn Zst /sV respectively. Consequently, T s „ st/st , is not a true coordinate 
transformation matrix in the rigorous sense. 
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Azimuth and coning angle misalignment effects must also be applied to the accelerometer 
sensor package associated with the model fixed reference frame. The development is highly 
similar to Equations (III-B.7)-(III-B.14) in that each axis must be treated individually, and a 
sequential transformation of first azimuth rotation and second coning rotation is considered. 
Details of this process will not be discussed but are expressed in Equations (III-B. 1 5)-(III-B.22). 
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Finally consider bias and sensitivity characteristics of the sensor packages. The sting fixed 
and model fixed accelerometer sensor package bias and sensitivity characteristics are modeled by 
the following relationships. 
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In Equations (III-B.23)-(III-B.24), V sn . denotes the sensor package output voltage vector, V x - 

i snj 

V v -V z denote components of V sn ., B sn . denotes the sensor package bias vector, B x -B v - 
y sn i ^snj i snj J sn| 

B Zsn denote components of B sn ., S sn . denotes the sensor package diagonal sensitivity matrix, 
S x -S y -S z denote components of S sn „ and i = st, m for the sting fixed and model fixed 
sensor packages. V s „ st and V s „ m are considered to be known vectors outputed by the sensors. 

Overall vector expressions for the accelerometer sensor packages are obtained by 
combining Equations (III-B.6), (III-B.14), and (III-B.23) and (III-B.5), (III-B.22), and (III-B.24). 
These expressions are 
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In scalar form, Equations (III-B.25)-(III-B.26) become 
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(III-B.32) 


where the notation C a and S a denote cos(a) and sin(a). Equations (III-B.27)-(III-B.32) provide 
the underlying analytical framework between the multi-axis gravitational acceleration 
measurement data in voltage form and the desired three dimensional Eulerian angular state of the 
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test model. These expressions contain both calibration parameters and the rotation angles, and 
provide the basis for accelerometer sensor calibration and attitude construction logic presented in 
Section IV. 


C. Gyro Equations 

Derivation of the underlying analytical framework between multi-axis angular rate data 
measured by a single rate gyro sensor package and the desired three dimensional Eulerian 
angular state of the test model is considered in this subsection. This analytical framework 
provides the basis for gyro sensor calibration and attitude construction logic presented in Section 
IV. With the geometry described in Figure 4, the angular velocity vector (algebraic form) of the 
model expressed in the model fixed reference frame is 




^ Oy ^ 

A m 7m z 


IT 
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(III-C.l) 


where Q Xm -Qy m ~Q Zm denote the model angular velocity components. Q m is considered to be a 
known vector outputed by the sensors (or determined from the sensor outputs). Extraction of the 
complete pitch-yaw-roll orientation state from an angular rate measurement and integration 
scheme employing a single gyro sensor package is feasible. Therefore, in the general framework 
of this research, angular velocity will only be measured in the model fixed axes using a single 
gyro sensor package. By relating these measured angular rate components with the temporal 
rates of change of the Euler rotation angles, the full angular state can be determined. This 
approach is applied to development of the general angular measurement system. A simplified 
case corresponding to the 31 -Inch Hypersonic Tunnel Application (fixed known roll angle) can 
always be extracted from the general framework. 

Based on the pitch-yaw-roll rotation sequence, the angular velocity vector can be built-up 
from the component rotation rates, or 

= ♦Tm/A + 'lm/A (III-C.2) 
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(I1I-C.3) 
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(III-C.5) 


Note Euler pitch rate 0 is executed about the first intermediate y axis (y j|), Euler yaw rate 4* is 
executed about the sting z axis (z st ), and Euler roll rate 4> is executed about the model x axis 
(x m ). In Equation (III-C.2), T,,,/,,, is the identity matrix, and Tm/jj is computed from the 
component transformation matrices T st /jj and 

Now consider effects from misalignment errors due to the azimuth and coning angles in the 

model mounted gyro sensor package. The following analysis must be conducted individually for 

each desired measurement axis. First consider the model x m axis as the desired measurement 

axis. A sequential transformation consisting first of azimuth rotation (A Xm ) about x m leading to a 

second intermediate frame, and next of coning rotation (C Xm ) about the intermediate y axis 

which leads to the true measurement or sensor fixed reference frame (x sn y sn z sn ) 

x m x m x m 

associated with x m , is implemented. For the y m axis, azimuth rotation (Ay m ) is about y m and 

coning rotation (C y ) is about the intermediate z axis leading to the y m specific sensor fixed axes 

(x sn ysn z sn )• Finally, for the z m axis, azimuth rotation (A z ) is about z m and coning 
ym >m ym m 

rotation (C z ) is about the intermediate x axis leading to the z m specific sensor fixed axes 
( x sn ysn z sn )• These transformation relationships are summarized in Equations (III-C.6)- 

z m z m z m 

(III-C. 11) where £2; is the model angular velocity expressed in frame i, Q x .-Q y .-Q z . are the 
components of £2j, and Tj/j is the transformation matrix from frame j to frame i. Here i, j = m, 
i2 k , sn k representing model fixed, second intermediate, and sensor fixed axes, respectively, and k 

= x m> y m . and z m- 

^i2 v = T i2 v /nAn 


‘m 


m 


(I1I-C.6) 
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L sn v M x 
m m 


cos(C x )0-sin(C ) 
x m x m 

0 1 0 

sin(C x ) 0 cos(C ) 
x m x m 


Q i2 v “ T i2 v /m^n 
•m ’m 


(III-C.8) 
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z m An 
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0 cos(C 7 ) sin(C 7 ) 
z m z m 


0-sin(C 7 )cos(C ) 
z m An 


The angular velocity vector expressed in 
the model fixed angular rate vector is obtained 



T 

sn v /m 

x m 


three sensor fixed reference frames in terms of 
combining Equations (II1-C.6)-(III-C. 1 1), or 
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(III-C.12) 


Q, 


sn>i 


= T. 


m 



An 


T 

sn 7 An 

L m 

The overall transformation matrices T sn /m , T sn /m , and T sn /m from model fixed to sensor 

x m -'m z m 

fixed axes in Equation (III-C.12) are obtained from a product of the component transformation 
matrices. Note only the x, y, and z components, respectively, in these three vector expressions 
are of interest. Extracting these three components and constructing a new vector relationship 
yields Equation (III-C.13). 


^n m ^sn m /m^m 




m 



T 

sn m An 


sn x m /m l 

rsn y m /m 2 

r 

l sn Zm An 3 


(III-C.13) 


In Equation (III-C.13), T s „ x / mi , T sny / m 2 , and T s „ z / m3 denote the first, second, and third 
rows of T sn /m , T sn /m , and T sn /m , respectively. Consequently. T sn /m , is not a true 
coordinate transformation matrix in the rigorous sense. 

Finally consider bias and sensitivity characteristics of the sensor package. The model fixed 
gyro sensor package bias and sensitivity characteristics are modeled by the following 
relationships. 


V sn - ®sn + ^sn ^sn 
sn m sn m sn m sn m 


r sn 


m 


B, 


sn 


m 


V V V 

Xsn m ysn m Zsn m 


B x B B 

sn m ysn m Zsn m 


'sn 
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0 


'sn 

0 

0 


m 


s v 

>sn. 


m 

0 s. 


0 

0 

z sn 


m 


(III-C.14) 


In Equation (III-C. 14), V sn denotes the sensor package output voltage vector, V x -V v 

m sn m ^sn m 

V z denote components of V sn , B sn denotes the sensor package bias vector, B x -B v 
sn m hi m sn m * sn jyj 

B Zcn denote components of B s „ m , S si , m denotes the sensor package diagonal sensitivity matrix, 
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and S x -S v -S 7 denote components of S«. n . V,. n is considered to be a known vector 
x sn m >sn m z sn m r s>n m sn m 

outputed by the sensor. 

An overall vector expression for the gyro sensor package is obtained by combining 
Equations (III-C.2), (III-C. 1 3), and (III-C.14). This expression is 

Ysn m = B sn m + ^snnjTsnjnAn^myil^il + ^m/st^st + ) (III-C. 15) 


In scalar form, Equation (III-C. 15) becomes 

V Xcn =B x en + S x sn ^ C C Y S V + S Av S C V C W C <D +C A V S Cv C W S <&)® 


' sn m Asn m ~ sn m v x m 


x m v ' x m 


x m x m 


+ (S Ax s c S*-C A s c C*)*+(C c w 

m m m m m 

v y sn =B y sn +s y S n ^ (_c Av s c v S w + C c v c t c <d~ s a v s c v c w s <t>) 0 


m 


m 


'y m y m 


7 m J m 


L sn 


+ (C C S (D + S A V S c v + (~ C A y S Cy 

= B z cn + S z sn « S A 7 S C 7 S »P _C A 7 S C 7 C^-Cc, Cq,S^)0 
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m 


m 


m m 


z m ^ 


m 


(III-C. 16) 


(III-C. 17) 


(III-C. 18) 


+ (-Ca Sp S^+Cp C*)#+(S* Sp )&} 


Z 7. 

m A m 


m 




- m 


Equations (III-C. 16)-(III-C. 18) provide the underlying analytical framework between the multi- 
axis angular rate measurement data in voltage form and the temporal rate of change of the 
desired three dimensional Eulerian angular state of the test model. These expressions contain 
both calibration parameters and the rotation angles and their time derivatives, and provide the 
basis for gyro sensor calibration and attitude construction logic presented in Section IV. 
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Section IV 

Sensor Calibration and Attitude Construction 


A. Accelerometer Sensor Calibration 

Before the AMS prototype can be used in the attitude construction mode, the accelerometer 
sensor packages must be calibrated. A consistent set of sensor output voltages and pitch-yaw-roll 
excitation angles, obtained from a controlled high-precision laboratory environment, can be used 
to determine unknown calibration parameters for each of the three accelerometers in both the 
sting fixed and model fixed packages. Synthetically generated data simulating the laboratory 
environment, however, will be used for calibration procedure development. Calibration 


parameters to be addressed include bias (B x , B v , B z , B x , B v , B z ), sensitivity 

sn st ^ sn st sn st sn m ' sn m sn m 

(S x , S v , S 7 , S x , S v , S 7 ), azimuth angle (A x . A v . A 7 . A x , A v , A 7 ), and 
coning angle (C Xst , Cy st , C Zst , C Xm , Cy m , C Zm ). The calibration parameters for each 
accelerometer are determined independently but in a sequential fashion. Sting fixed 


accelerometer sensors are addressed first, and model fixed accelerometer sensors are considered 


second. If deemed to be important, other sensor imperfections could be incorporated into the 
calibration procedure. For example, temperature corrections, although judged to be an important 
factor for the 31 -Inch Hypersonic Tunnel application, were not considered in the prototype 
development stage. However, calibration for thermal effects can be appended to the algorithm 
presented in this section at a more advanced AMS development stage. 

A Gauss-Newton iterative nonlinear regression method outlined in Reference 5 will be 


used as a basis for accelerometer calibration. Further details of the Gauss-Newton algorithm can 
be found in Reference 17. Recall the governing sting fixed accelerometer sensor relationships in 
Equations (III-B.27)-(III-B.29). From a sensor calibration perspective, voltages, angular 
orientation states, and gravity are known parameters while bias, sensitivity, azimuth angle, and 


coning angle are unknown variables implying a set of nonlinear algebraic equations must be 
solved. Implications from sensor model assumptions are that each scalar equation involves only a 
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single unknown quadruple (bias, sensitivity, azimuth angle, coning angle), allowing an 
independent solution strategy. Represent any one of the three sting fixed accelerometer sensor 
relationships as 

V=F(B,S,A,C,0,T\G) (IV-A.l) 


where function F denotes the nonlinear functional dependency of voltage V on bias B, sensitivity 
S, azimuth rotation A, coning rotation C, pitch-yaw angles ©-'F, and gravity G. The right-hand 
side of this equality can also be expanded in a Taylor series about a reference solution. 

v = v r + dF fx P \ (X - X r ) + ... (IV-A.2) 


x = [b s a c] t 


dF(X,P) _[ dF BF BF BF~ 
c)X c)B d A dC 


p = [e>p g] t (IV . a . 3) 


In Equation (IV-A.2), X denotes a vector containing the unknown variables, P denotes a vector 
of known parameters, and subscript "r" denotes evaluation at the known reference condition. 

Suppose initial estimates for the bias, sensitivity, azimuth, and coning variables are known 
and interpreted to be the reference solution X r in Equation (IV-A.2). Further suppose that a 
consistent data set of n quadruples for voltage, pitch, yaw, and gravity is available. A single 
quadruple from this set can be interpreted as V and P appearing in Equation (IV-A.2). Using the 
reference solution and parameter data, voltage and it’s partial derivatives can also be computed 
from Equation (IV-A.l) and interpreted as V r and dF(X,P)/dX r . Construct a vector equivalent to 
Equation (IV-A.2) using the nx( 1+2+1) data set to yield 

V = v r + aF ^ — (X -X r ) + ... (IV-A.4) 
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5 B 2r 


dF dF 


dS ] dA ] 


dF dF 
dS? dA 7 


OF 

dC ] 
*r 

dF 

dC t 


OF dF dF dF 
^n,. ^n r ^'n r ^ r n r 


(IV-A.5) 


Multiply this expression by 6F(X,P)/t)X r T on the left, perform matrix inversion, neglect higher 
order terms, and rearrange to yield 


. ,dF(X,P)T 5 F(X,P) ,_,aF(X,P)T^ 

X-X r +( 5X r 3X r } ax r (V_V r ) 


(IV-A.6) 


Equation (IV-A.6) represents the iterative Gauss-Newton regression algorithm. X is computed 
from Equation (IV-A.6) and reinterpreted as the reference solution for the next computational 
cycle. This process is repeated until a solution tolerance is met. 

For completeness, partial derivatives of Equations (III-B.27)-(III-B.29) are given below. 
For the sting x stlst axis accelerometer sensor, partial derivatives are 



(IV-A.7) 
(IV-A.8) 
(IV-A.9) 
(IV- A. 10) 


For the sting y sr , st axis accelerometer sensor, partial derivatives are 
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(IV -A. 11) 
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For the sting z sns( axis accelerometer sensor, partial derivatives are 
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(IV-A.12) 

(IV-A.13) 

(IV-A.14) 

(IV-A.15) 

(IV-A.16) 

(IV-A.17) 

(IV-A.18) 


Now focus attention on calibration of the model fixed accelerometer sensors and their 
corresponding relationships in Equations (III-B.30)-(II1-B.32). Development is highly similar to 
Equations (IV-A.1)-(IV-A.18) and is only briefly discussed here. Voltages, angular orientation 
states (d> only), and now sting fixed gravity components are the known parameters while bias, 
sensitivity, azimuth angle, and coning angle remain the unknown variables leading to a set of 
nonlinear algebraic equations. Note the sting fixed gravity components can be computed from 
Equation (III-B.6). The iterative Gauss-Newton regression algorithm represented in Equations 
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(IV-A.1)-(IV-A.6) is again applicable with the only changes affecting the functional dependence 
on parameters. For the model fixed sensor package, expressions for V and P become 


V = F(B, S,A,C,*,G Xst ,G yrf G ZBi ) 


P = 


<FG X G v G 7 

x st y s t z st 


it 


(IV-A.19) 

(1V-A.20) 


For completeness, partial derivatives of Equations (III-B.30)-(III-B.32) are given below. 
For the model x snm axis accelerometer sensor, partial derivatives are 
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(IV-A.21) 


(1V-A.22) 
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(IV- A. 24) 


For the model y sn axis accelerometer sensor, partial derivatives are 
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(IV-A.27) 
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For the model z sn axis accelerometer sensor, partial derivatives are 
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(IV-A.32) 


The Gauss-Newton regression algorithms for accelerometer sensor calibration are validated 
here within a simulated laboratory environment with synthetically generated data using 
MATLAB. It is important to note the general calibration algorithms presented here are not tested 
in the laboratory with real data within the data acquisition control software environment using 
LabVIEW. An angular motion profile was generated and the corresponding sensor voltages were 
computed for specified bias, sensitivity, azimuth angle, and coning angle parameters. This data is 
intended to represent actual data that would be collected from hardware using the angular 
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calibration table in the laboratory. The motion profile consisted of a linear progression from 0 
deg to 10 deg for 0- l f , -d> over 10 s with a time step of 0.01 s which yields n = 1,001. Table 3 
lists the specified and computed calibration parameters for all accelerometer sensors. Specified 
values are based on the manufacture’s published data. The computed values are all in excellent 
agreement with the specified values validating the accelerometer sensor calibration 
methodology. Appendix B lists the MATLAB software pertaining to accelerometer sensor 
calibration. 


Table 3 Accelerometer Calibration Test Results 


Axis 

Parameter 

Specified 

Computed 

Specified 

Computed 


B (v) 

S (v/g) 

x sn st 

1.25 

1.2500 

0.51 

0.5100 

y S n st 

1.27 

1.2700 

0.47 

0.4700 

z sn st 

1.19 

1.1900 

0.54 

0.5400 

Xsn m 

1.26 

1.2600 

0.52 

0.5200 

^ sn m 

1.28 

1.2800 

0.49 

0.4900 

Zsn m 

1.20 

1.2000 

0.53 

0.5300 


A (( 

leg) 

C (c 

•eg) 

x sn st 

0.05 

0.0500 

0.06 

0.0600 

ysn st 

0.12 

0.1200 

0.15 

0.1500 

! 

z sn st 

0.11 

0.1100 

0.11 

0.1100 

x sn m 

0.08 

0.0800 

0.08 

0.0800 

! 

y S n m 

0.10 

0.1000 

0.09 

0.0900 

z snm 

0.09 

i 

1 

0.0900 

0.12 

0.1200 
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B. Accelerometer Attitude Construction 

Once the accelerometer sensor packages have been calibrated, the AMS prototype can be 
used for attitude construction. Given numeric knowledge of the calibration parameters, a set of 
sensor output voltages can be used to determine unknown pitch-yaw-roll angles associated with 
the two accelerometer packages mounted in the sting fixed and model fixed reference frames. To 
verify accurate attitude construction, utilization of a consistent set of sensor output voltages and 
pitch-yaw-roll excitation angles, obtained from a controlled high-precision laboratory 
environment, is highly desirable. Synthetically generated data simulating the laboratory 
environment, however, will be used for construction procedure development. The attitude states 
are constructed sequentially. Pitch-yaw angles (0-W) are addressed first, and roll angle (d>) is 
considered second. Two different construction methodologies are considered and include an 
analytic closed-form strategy and an iterative numeric strategy. The closed-form strategy 
facilitates a real-time attitude construction submode for the AMS, while the iterative strategy is 
incrementally more accurate. 5 

The analytic closed-form attitude construction strategy will be addressed first. Recall the 
governing sting fixed and model fixed accelerometer sensor relationships in Equations (III- 
B.25)-(1II-B.26). After performing matrix operations on these relationships, one can find 

T sn st /st^sn st ^sn st " B sn st > = T st/t G t (IV-B. 1 ) 

^sn m /m^sn m ^sn m - ®sn m ) = ^m/st G st (IV-B. 2) 

Using Equations (II1-B.5)-(I1I-B.6), note the left-hand side of Equation (IV-B.l) can be renamed 
as G st , while the left-hand side of Equation (IV-B. 2) can be replaced with G m . In scalar form. 
Equations (IV-B. 1)-(IV-B.2) become 



°Z S , = c e° 

(IV-B.5) 

G = G, 

x m x st 

(IV-B.6) 

j v = C^G -f-S^G 
ym <p yst ® z st 

(IV-B.7) 

z m ° yst •P z st 

(IV-B.8) 


From an attitude construction perspective, all gravity terms are known parameters while pitch- 
yaw-roll angles are unknown variables implying a set of nonlinear algebraic equations must be 
solved. 

Manipulation of Equations (IV-B.3)-(IV-B.5) to isolate © and W yields 


JgI +Gj G x 

© =tan -l (— — g- is*-) S gn ( -- ^4- -£*) , 'F * ±90 deg 


^St 


W = tan" 1 (- G ^) 

x st 


(IV-B.9) 

(IV-B.10) 


Further manipulation with Equations (IV-B.7)-(IV-B.8) to isolate O yields 

ym y s t z m z st 

The sign function in Equation (IV-B.9) is appended to the solution to recover directional pitch 
information. Given the gravity components. Equations (IV-B.9)-(IV-B.l 1) offer analytic closed- 
form solutions for attitude construction. It is important to note that gravity components are 
computed from the sensor output voltages and calibration parameters appearing on the left-hand 
sides of Equations (1V-B.1)-(1V-B.2). The analytic closed-form solution strategy is highly 
attractive in applications requiring real-time monitoring of wind tunnel testing. 

Now consider the iterative numeric attitude construction strategy. The Gauss-Newton 
regression method will once again be applied here as outlined in Reference 5. Recall the 
governing sting fixed accelerometer sensor relationships in Equations (I1I-B.27)-(III-B.29). From 
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an attitude construction perspective, voltages, calibration terms, and gravity are known 
parameters while pitch-yaw angles are unknown variables implying a set of nonlinear algebraic 
equations must be solved. Mathematical pitch-yaw coupling structure implies all three 
relationships must be considered simultaneously in the solution strategy since each scalar 
equation involves the same unknown pair (pitch angle, yaw angle). Represent any one of the 
three sting fixed accelerometer sensor relationships as 

Vj = Fj(0 , V, Bj, Sj, Aj, q, G) (IV-B.12) 


where function Fj denotes the nonlinear functional dependency of voltage Vj on pitch-yaw angles 
0-W, bias Bj, sensitivity Sj, azimuth rotation Aj, coning rotation C,, and gravity G where i = 
x SIlst , y S n s p z sn st - The right-hand side of this equality can also be expanded in a Taylor series 
about a reference solution. 

Vj =Vj r + aF jx Pj) (X-X r ) + ... (IV-B.13) 


X = [0 yf 


3Fj(X,Pj) 

dx 


3E 3Fj 
30 W 


Pj= Bj Sj AjCj g] T (IV-B.14) 


In Equation (IV-B.13), X denotes a vector containing the unknown variables, P; denotes a vector 
of known parameters, and subscript "r" denotes evaluation at the known reference condition. 

Suppose initial estimates for the pitch-yaw variables are known and interpreted to be the 
reference solution X r in Equation (IV-B.13). Further suppose that a data set of n triples for each 
voltage is available. A single triple from this set can be interpreted as Vj appearing in Equation 
(IV-B.13) for i = x stlst , y snst , z sns( . Further, P; is known from calibration processing and the 
environment. Using the reference solution and parameter data, voltages and their partial 
derivatives can also be computed from Equation (IV-B.12) and interpreted as V if and 
dFj(X,P,)/dX r . Construct a vector equivalent to Equation (IV-B.13) using the nx3( 1+4+1) data 
set to yield 
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(IV-B.15) 
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(IV-B.16) 


Multiply this expression by 6F(X,P)/dX r T on the left, perform matrix inversion, neglect higher 
order terms, and rearrange to yield 


, ,dF(x,P)TaF(x,P) N _!aF(x,P)T /ltr ^ 

x_x r +(- 5X ■ r 5X ) ax r (V ~v 


(IV-B.17) 


Equation (1V-B.17) represents the iterative Gauss-Newton regression algorithm. X is computed 
from Equation (1V-B.17) and reinterpreted as the reference solution for the next computational 
cycle. This process is repeated until a solution tolerance is met. 

For completeness, partial derivatives of Equations (III-B.27)-(III-B.29) are given below. 
For the sting x s „ st axis accelerometer sensor, partial derivatives are 
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av x 

A $ jq 

7777 = ^ x {Gp SqSuz+Sa S c SqCu;}G 
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(IV-B.19) 


For the sting y sr , st axis accelerometer sensor, partial derivatives are 
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(IV-B.20) 
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For the sting z sr , st axis accelerometer sensor, partial derivatives are 

dV, 


L sn 
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Now focus attention on attitude construction using the model fixed accelerometer sensors 
and their corresponding relationships in Equations (III-B.30)-(III-B.32). Development is highly 
similar to Equations (IV-B. 12)-(IV-B.23) and is only briefly discussed here. Voltages, 
calibration terms, and now sting fixed gravity components are the known parameters while roll 
angle is the unknown variable leading to a set of nonlinear algebraic equations. Note the sting 
fixed gravity components can be computed from Equation (III-B.6). The iterative Gauss-Newton 
regression algorithm represented by Equations (IV-B.12)-(IV-B.17) is again applicable with the 
only changes affecting the unknown variable and known parameter definitions, and the 
functional dependence on these quantities. For the model fixed sensor package, expressions for 
Vj, X and Pj become 


Vj Fj(d>, Bj, Sj, A j, Cj, G x ^, Gy^, G z ^) 


(IV-B. 24) 
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For completeness, partial derivatives of Equations (III-B.30)-(III-B.32) are given below. 
For the model x sr , m axis accelerometer sensor, partial derivatives are 


44 



(IV-B.26) 
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For the model y s „ m axis accelerometer sensor, partial derivatives are 
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For the model z sr , m axis accelerometer sensor, partial derivatives are 
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(IV-B.28) 


The analytic closed-form and iterative numeric algorithms for accelerometer attitude 
construction are validated here within a simulated laboratory environment with synthetically 
generated data using MATLAB. It is important to note the general construction algorithms 
presented here are not tested in the laboratory with real data within the data acquisition control 
software environment using LabVIEW. An angular motion profile was generated and the 
corresponding sensor voltages were computed for specified bias, sensitivity, azimuth angle, and 
coning angle parameters. This data is intended to represent actual data that would be collected 
from hardware using the angular calibration table in the laboratory. The motion profile consisted 
of a linear progression from 0 deg to 10 deg for 0-W-<I> over 10 s with a time step of 0.01 s 
which yields n = 1,001. Specified calibration parameter values from Section IV-A are used. 
Figure 6 shows overlay traces of the exact and computed angular states (both analytic closed- 
form and iterative numeric strategies). Figures 7-8 show difference traces between exact and 
computed angular states. The computed values are all in excellent agreement with the exact 
values validating the accelerometer sensor calibration methodology. Figures 7-8 suggest the 
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analytic closed-form solution is slightly more accurate than the iterative numeric solution. With 
noise and resolution effects considered, the iterative numeric solution would be expected slightly 
more accurate. Appendix C lists the MATLAB software pertaining to accelerometer attitude 
construction. 
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Figure 6 Accelerometer Construction Test Results - Absolute 


47 





IPhi Errorl (deg) IPsi Error! (deg) ITheta Error! (deg) 



i i i i i i i i i 1 1 

0 1 2 3 4 5 6 7 8 910 


t (s) 


Figure 7 Accelerometer Construction Test Results - Relative (Analytic) 
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C. Gyro Sensor Calibration 

Before the AMS prototype can be used in the attitude construction mode, the gyro sensor 
package must be calibrated. A consistent set of sensor output voltages and pitch-yaw-roll 
excitation angles and angle rates, obtained from a controlled high-precision laboratory 
environment, can be used to determine unknown calibration parameters for each of the three 
gyros in the model fixed package. Synthetically generated data simulating the laboratory 
environment, however, will be used for calibration procedure development. Calibration 
parameters to be addressed include bias (B x , B v , B z ), sensitivity (S x , S v , S z ), 
azimuth angle (A Xm , Ay m , A Zm ), and coning angle (C Xm , Cy m , C Zm ). The calibration parameters 
for each gyro are determined independently. If deemed to be important, other sensor 
imperfections could be incorporated into the calibration procedure. For example, temperature 
corrections, although judged to be an important factor for the 31-Inch Hypersonic Tunnel 
application, were not considered in the prototype development stage. However, calibration for 
thermal effects can be appended to the algorithm presented in this section at a more advanced 
AMS development stage. 

A Gauss-Newton iterative nonlinear regression method will be used as a basis for gyro 
calibration. Further details of the Gauss-Newton algorithm can be found in Reference 17. Recall 
the governing model fixed gyro sensor relationships in Equations (III-C.16)-(III-C. 18). From a 
sensor calibration perspective, voltages, angular orientation states, and angular orientation state 
derivatives are known parameters while bias, sensitivity, azimuth angle, and coning angle are 
unknown variables implying a set of nonlinear algebraic equations must be solved. Implications 
from sensor model assumptions are that each scalar equation involves only a single unknown 
quadruple (bias, sensitivity, azimuth angle, coning angle), allowing an independent solution 
strategy. Represent any one of the three model fixed gyro sensor relationships as 


V =F(B, S, A,C, *P, d>, 0, <I>) 


(IV-C.l) 
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where function F denotes the nonlinear functional dependency of voltage V on bias B, sensitivity 
S, azimuth rotation A, coning rotation C, yaw-roll angles 'F-O, and pitch-yaw-roll angle rates 
@_vp_<j> The right-hand side of this equality can be expanded in a Taylor series about a 
reference solution. 


v = V r + - F( 3 x P) r (X-X r ) + . 


x = [b s a c] T = 


dF dF dF dF 

aB ^ aA 


p = 


vp<d e'i'fi) 11 


(IV-C.2) 

(IV-C.3) 


In Equation (IV-C.2), X denotes a vector containing the unknown variables, P denotes a vector 
of known parameters, and subscript "r" denotes evaluation at the known reference condition. 

Suppose initial estimates for the bias, sensitivity, azimuth, and coning variables are known 
and interpreted to be the reference solution X r in Equation (IV-C.2). Further suppose that a 
consistent data set of n sextuples for voltage, yaw, roll, pitch rate, yaw rate, and roll rate is 
available. A single sextuple from this set can be interpreted as V and P appearing in Equation 
(IV-C.2). Using the reference solution and parameter data, voltage and it’s partial derivatives can 
also be computed from Equation (IV-C.l) and interpreted as V r and dF(X,P)/dX r . Construct a 
vector equivalent to Equation (IV-C.2) using the nx( 1+2+3) data set to yield 


V=V r + 
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(IV-C.5) 
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Multiply this expression by dF(X,P)/6X r T on the left, perform matrix inversion, neglect higher 
order terms, and rearrange to yield 


v _ v . ,dF(XJ)TdF(X,P) ,_,aF(X4»)T^ y 
X_X r +( 3X- r - ax J - 5X~ r (V-V r> 


(1V-C.6) 


Equation (IV-C.6) represents the iterative Gauss-Newton regression algorithm. X is computed 
from Equation (IV-C.6) and reinterpreted as the reference solution for the next computational 
cycle. This process is repeated until a solution tolerance is met. 

For completeness, partial derivatives of Equations (III-C. 16)-(III-C.18) are given below. 
For the model x SI , m axis gyro sensor, partial derivatives are 
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For the model y srim axis gyro sensor, partial derivatives are 
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(IV-C.12) 
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For the model z SIlm axis gyro sensor, partial derivatives are 
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The Gauss-Newton regression algorithm for gyro sensor calibration is validated here within 
a simulated laboratory environment with synthetically generated data using MATLAB. It is 
important to note the general calibration algorithm presented here is not tested in the laboratory 
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with real data within the data acquisition control software environment using LabVIEW. An 
angular motion profile was generated and the corresponding sensor voltages were computed for 
specified bias, sensitivity, azimuth angle, and coning angle parameters. This data is intended to 
represent actual data that would be collected from hardware using a high bandwidth angular 
calibration table in the laboratory. The motion profile consisted of harmonic excitation with 0 = 
10sin(0.125x2jrxt), W = 10sin(0.0625x2jtxt), and O = 10sin(0.25x2jtxt) deg for 0 s t < 10 s with 
a time step of 0.01 s which yields n = 1,001. Table 4 lists the specified and computed calibration 
parameters for all gyro sensors. The computed values are all in excellent agreement with the 
specified values validating the gyro sensor calibration methodology. Appendix D lists the 
MATLAB software pertaining to gyro sensor calibration. 


Table 4 Gyro Calibration Test Results 


Axis 

Parameter 

Specified 

Computed 

Specified 

Computed 


B (v) 

S (v/c 

eg/s) 

x sn m 

0.051 

0.05100 

1.25 

1.2500 

ysn m 

0.052 

0.05200 

1.26 

1.2600 

Zsn m 

0.053 

0.05300 

1.27 

1.2700 


A (c 

leg) 

C (c 

leg) 

Xsn m 

0.14 

0.1400 

o.ii 

0.1100 

ysn m 

0.15 

0.1500 

0.12 

0.1200 

Zsn m 

0.16 

0.1600 

0.13 

0.1300 


D. Gyro Attitude Construction 

Once the gyro sensor package has been calibrated, the AMS prototype can be used for 
attitude construction. Given numeric knowledge of the calibration parameters, a set of sensor 
output voltages can be used to determine unknown pitch-yaw-roll angles associated with the 
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single gyro package mounted in the model fixed reference frame. To verify accurate attitude 
construction, utilization of a consistent set of sensor output voltages and pitch-yaw-roll 
excitation angles and angle rates, obtained from a controlled high-precision laboratory 
environment, is highly desirable. Synthetically generated data simulating the laboratory 
environment, however, will be used for construction procedure development. All attitude states 
(G-W-O) are constructed in a simultaneous fashion due to the coupled gyroscopic dynamics 
associated with this class of sensors. Analytic closed-form solutions are not feasible here, thus 
leaving only iterative numeric strategies for gyro attitude construction. 

Recall the governing model fixed gyro sensor relationship in Equation (III-C.15). After 
performing matrix operations on this relationship, one can find 

^sn m /m^sn m ^sn m - ®sn m ^ = ^m/il® il + ^m/st^st + ^m/m^m (IV-D. 1 ) 


Using Equations (III-C.3)-(III-C.5), note the right-hand side of Equation (IV-D.l) can be 
compressed into a product of a single matrix and vector, or 


<A a = [© 'f' <t>] 


m/il®il + ^m/st^st + ^ ml eft a 

|T 


T = 
m/a 


T T T 
m/il2 m/ist3 m/m j 


(IV-D.2) 

(IV-D.3) 


In Equation (IV-D.3), Tm/stj, and T,n/ m j denote the first, second, and third columns of 

Tm/n> T jn/gj, and T respectively. Consequently, T,,^, is not a true coordinate 
transformation matrix in the rigorous sense. Notational usage of the letters A and a denote angle. 
Combining Equations (IV-D.1)-(IV-D.2) with further matrix manipulation yields 

a «=' I 'm 1 ,a T „ m /„,S^„/V s „ in -B sniiil (IV-D.4) 

Examination of Equation (IV-D.4) reveals that the right-hand side is a function of the angular 
states, sensor output voltages, and sensor calibration parameters, while the left-hand side is a 
function of the angular state derivatives. From an attitude construction perspective, voltages and 
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calibration terms are known parameters while pitch-yaw-roll angles and their derivatives are 
unknown variables implying a set of nonlinear coupled differential equations must be solved. 
Gyro attitude construction is thus a fundamentally different problem than accelerometer based 
attitude construction. 

A 5th order Runge-Kutta numerical integration method will be used as a basis for gyro 
construction. Details of the Runge-Kutta algorithm can be found in Reference 17. This algorithm 
requires differential equations in 1st order form. Note Equation (IV-D.4) satisfies this 

i 

requirement. Represent any one of the three model fixed gyro sensor relationships within 
Equation (IV-D.4) as 

a i = FiCA j , V k , B k , S k , A k , C k ) < IV - D ' 5 > 

I 

where function Fj denotes the nonlinear functional dependency of angular state derivative 9lj on 

! 

angular state 9lj, voltage V k , bias B k , sensitivity S k , azimuth rotation A k , and coning rotation C k , 
where i,j = 1, 2, 3, 91] = 0, 9l 2 = W, 9l 3 = 3>, and k = x sn , y S n m ’ z sn m - Equation (IV-D.5) can be 
interpreted as a nonlinear dynamic system with 9lj as state variables and V k as inputs. 

Suppose initial conditions for the pitch-yaw-roll angular states are known from the pretest 
| AMS alignment. Further suppose that a data set of n triples for each voltage is available. With 

I 

the input voltages, angular position starting conditions, and calibration parameters, Equation (IV- 
D.5) can be used to generate state derivatives needed for temporal propagation of the angular 
states through numerical integration. In the 5th order Runge-Kutta strategy, denote the current 
angular state at time t as 9lj(t) and let the time step be denoted by At. To advance the state to 
5lj(t+At), Equation (IV-D.5) is evaluated 6 times at sequentially projected preliminary angular 
state solutions across the time step, or 
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(IV-D.6) 


j - Fj(^j(t), V k (t), B k , S k , A k , C k ) 
ft j 2 = Fj(9l j (t) + j j At, V k (t + ^At), B k , S k , A k , C k ) 

= Fj(5lj(t) + j] +! ^ j2^ t ’ ^k(t + ^At), B k , S k , A k , C k ) 

^i 4 = Fj(9lj(t) + ^-{-91 j 2 +29l j 3 )At, V k (t + ^-At), B k , S k , A k , C k ) 

&i 5 = Fj(9lj(t) + -^-{$1 ji+3$l j 4 )At, V k (t + ^-At), B k , S k , A k , C k ) 

A i6 = F i (?l j (t) + i{-3!Aj 1 +2!A j2 +12!A j3 -12!?l j4 +8^lj 5 }At, V k (t + At), B k , S k , A k , C k ) 


The final estimate for the angular state is obtained by averaging 5 of the preliminary solutions 
according to 

9lj(t + At) = JljCt) + ^ J {7A il +32A i +1251^+32^5 +7?l i }At (IV-D.7) 

Equations (IV-D.6)-(IV-D.7) represent the iterative Runge-Kutta integration algorithm. ?lj(t+At) 
is computed from Equation (IV-D.7) and reinterpreted as the starting condition for the next 
computational cycle. This process is repeated until the final time is reached. The angular 
construction accuracy will be a direct function of both the time step size and the angular state 
initial condition. 

The numerical integration algorithm for gyro attitude construction is validated here within 
a simulated laboratory environment with synthetically generated data using MATLAB. It is 
important to note the general construction algorithm presented here is not tested in the laboratory 
with real data within the data acquisition control software environment using LabVIEW. An 
angular motion profile was generated and the corresponding sensor voltages were computed for 
specified bias, sensitivity, azimuth angle, and coning angle parameters. This data is intended to 
represent actual data that would be collected from hardware using a high bandwidth angular 
calibration table in the laboratory. The motion profile consisted of harmonic excitation with © = 
10sin(0.125x2jtxt), V = 10sin(0.0625x2:rtxt), and d> = 10sin(0.25x2:txt) deg for 0 < t < 10 s with 
a time step of 0.001 s which yields n = 10,001. Specified calibration parameter values from 
Section IV-C are used. Figure 9 shows overlay traces of the exact and computed angular states. 
Figure 10 shows difference traces between exact and computed angular states. The computed 
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values are all in excellent agreement with the exact values validating the gyro attitude 
construction methodology. Appendix E lists the MATLAB software pertaining to gyro attitude 
construction. 
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Figure 10 Gyro Construction Test Results - Relative 
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Section V 

Hardware-Software Design and System Integration-Test 


A. Hardware Description 

For the AMS prototype, very little freedom existed in hardware design or selection. No 
budgetary resources within the contract, or external to the contract, were available for hardware 
purchase such as MEMS sensors or data acquisition microprocessor cards. Therefore, utilization 
of existing NASA LaRC hardware inventories was mandated, and in particular inventories 
allocated to the Models and Systems Branch. Major hardware components that were made 
available to the contractors consisted of a Pentium II class personal computer with video 
monitor, National Instruments 4 channel 200 kHz NI-4452 data acquisition card. Analog Devices 
ADXL105 MEMS micromachined polysilicon vibrating structure accelerometer, and muRata 
ENC-03J piezoelectric vibrating gyroscope or rate gyro. Also available were several high fidelity 
angular calibration tables with motorized drive systems or manual cranks, instrumentation, 
cables and connectors, and laboratory space to conduct research activities. Critical sensor and 
data acquisition card components will be briefly discussed. 

An Analog Devices ADXL105 accelerometer was used for prototype AMS laboratory 
development activities based on availability considerations. Figure 1 1 shows a 3-view drawing 
of this MEMS accelerometer unit. The ADXL105 MEMS accelerometer uses an open-loop 
architecture consisting of a micromachined polysilicon vibrating structure outfitted with a central 
capacitance plate and positioned between two other fixed capacitor plates. The fixed plates are 
used to nominally vibrate the central plate, as well as to sense unsymmetric motion resulting 
from accelerations. The sensor outputs a voltage signal proportional to the input acceleration. 
Table 1 in Section II-C summarized the performance characteristics of this class of sensor, 
according to vendor published information. Detail sensor specifications can be found in 
Reference 13. 
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Figure 1 1 Analog Devices ADXL105 3-View Drawing (mm) 


A muRata ENC-03J gyroscope was also used for prototype AMS laboratory development 
activities based on availability considerations. Figure 12 shows a 3-view drawing of this MEMS 
gyro unit. The gyroscope is a rate gyro which produces an output voltage signal proportional to 
input angular rate. The ENC-03J MEMS gyro uses an open-loop architecture consisting of a 
ceramic bimorph vibrating structure excited electrically. When this oscillating unit experiences 
an angular rate, an electrical signal is generated proportional to the resulting Coriolis force. Table 
2 in Section II-C summarized the performance characteristics of this class of sensor, according to 
vendor published information. Detail sensor specifications can be found in Reference 15. 
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Figure 12 muRata ENC-03J 3-View Drawing (mm) 



A National Instruments NI-4452 data acquisition microprocessor card was available for 
prototype AMS laboratory development activity. Figure 13 shows a top-view photograph of this 
microprocessor card. The NI-4452 data acquisition card is a 16 bit resolution, 4 analog input 
channel, high throughput, high dynamic range signal acquisition device. Sampling rate for this 
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card is 200 kHz with an operational frequency range from DC to 95 kHz. The card incorporates 
sophisticated real-time anti-aliasing and anti-imaging digital signal processing filters. The NI- 
4452 also interfaces with lower level NI-DAQ driver software and higher level LabVIEW 
control software. Detail data acquisition card specifications can be found in Reference 18. 



Figure 13 National Instruments NI-4452 Top-View Photograph 

B. Software Design 

A wide range of design freedom and selection existed with data acquisition controlling 
software, such as C, Visual Basic, LabVIEW, MATLAB, Data Acquisition Toolbox, etc. 
National Instruments LabVIEW was selected for this function based on its powerful capabilities, 
ease of use, compatibility with the NI-4452 data card, and acceptance at NASA LaRC. 
LabVIEW’s graphical user interface was used to design multi-functional control instrumentation 
software that reads input signals from the data card, processes and filters the input signals, 
calibrates several key sensor parameters, transforms the acceleration and angular rate data to 
angular position states, records data to storage devices, and outputs real-time information to a 
video display in numeric and graphical formats. For the AMS prototype, and with a limited 
number of data acquisition channels imposed by the hardware, only simplified single-axis pitch 
calibration-construction software design using a single accelerometer and single gyro was 
considered in detail. A more capable multi-accelerometer based two angle construction 
capability was also designed, but to lesser degree of detail. Software design and testing was a 
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major effort within the overall project. Detail information about LabVIEW capabilities and 
utilization can be found in References 11-12. 

Figures 14-15 show the functional diagram of the AMS prototype control instrumentation 
software. In the LabVIEW environment, software design consists of graphically interconnecting 
system supplied and user defined component functions, in an efficient and effective arrangement, 
to accomplish higher level objectives, such as data acquisition or data processing. The AMS 
prototype operates in one of two modes (see Figure 3): sensor calibration or attitude 
construction, and Figures 14-15 reflect this architecture. Figure 14 pertains to sensor calibration 
while Figure 15 relates to attitude construction. These functional diagrams, and associated 
component level function diagrams, are the software listing, since LabVIEW auto-codes the low 
level commands. 

In the sensor calibration mode (see Figure 14), 3 analog signals are sampled by the 
software through the data acquisition card. These signals correspond to accelerometer output 
voltage, gyro output voltage, and a temperature voltage signal for the accelerometer (recall the 
ADXL105 unit has a built in temperature sensor). A significant problem was encountered during 
interfacing the NI-4452 card with core data sampling functions within LabVIEW. Success was 
ultimately achieved by utilizing and modifying several core data sampling functions from Wyle 
Laboratories. These component functions from Figure 14 are expanded in Figure 16. Next the 3 
element vector signal is split into separate scalar signals and run through low pass Chebyshev 
filters to reduce high frequency noise. Figure 17 shows the expanded functional diagram for the 
filters. The signals are then averaged over a sample window (assuming the sensors are at rest) to 
calibrate the bias parameters. In this preliminary development stage, calibration of sensitivity 
parameters is not considered. These bias values are then supplied to the attitude construction 
mode diagram in Figure 15 as the 3 signals entering from the right-hand side. 

In the attitude construction mode (see Figure 15), 3 analog signals are again sampled and 
correspond to accelerometer output voltage, gyro output voltage, and temperature voltage for the 
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Figure 14 Functional Diagram for Sensor Calibration 
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Figure 15 Functional Diagram for Attitude Construction 







accelerometer (assuming the sensors are in motion or have been moved from their initial state). 
The signals are sampled with the modified Wyle Laboratories core functions shown in Figure 16. 
As shown in Figure 15, the vector signal is again split and filtered in an identical fashion as that 
noted in Figure 14. For the accelerometer based attitude construction noted in Figure 15, scale 
factor (or sensitivity) and bias are applied to the accelerometer voltage signal. Note the 
sensitivity must be manually inputed by the user at this time. The resulting signal (gravity 
component measurement) is passed to an inverse sine function to generate pitch angle (see 
Equation (IV-B.3) with W = 0 deg), and then units are changed and the result is sent to the user 
interface panel for display. The gravity component measurement signal can also be passed to a 
more complex two angle attitude construction algorithm, as indicated in Figure 15. Figure 18 
shows the expanded functional diagram representing this algorithm. At this time, the user can 
manually input two additional gravity component signals and initial angle estimates to 
reconstruct the associated pitch-roll angles in a framework similar to that presented in Sections 
II1-B and IV-B. Resulting solutions are outputed to the user display panel. 

Gyro attitude construction logic noted in Figure 15 is structured differently. After 
subtracting off the bias value from the filtered gyro voltage signal, the result is temporarily 
stored in small packets and then numerically integrated as a unit with an initial user supplied 
value for the integration bin. After this processing, scale factor is taken into account and a default 
angle value is applied to generate a measured pitch angle, which is also displayed to the user. 
Temperature voltage signal is also processed with its appropriate calibration values resulting in a 
temperature measurement which could be used in scheduling accelerometer calibration 
parameters as the thermal environment changes. At this time, the temperature measurement 
signal is only displayed as an output. All control instrumentation and measurement functions 
displayed in Figures 14-15 can be performed in real-time allowing wind tunnel technicians and 
experimental aerodynamicists to evaluate ongoing test conditions. 
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Figure 16 Functional Diagram for Data Sampling 
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Figure 17 Functional Diagram for Data Filtering 
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Figure 18 Fiuictional Diagram for Accelerometer Two Angle Construction 
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Figure 19 User Interface Panel 


The user interface panel is shown in Figure 19. The top-left window in the interface panel 
displays the pitch angle measured by the accelerometer, and the top-right window displays the 
pitch angle measured by the gyro. Numeric displays labeled "Angle" in these two upper windows 
indicate measurement angles presented in real-time. Recent time histories of the pitch 
measurements are also displayed in graphical format. The bottom-left window also displays the 
two angle solution results. The mid-left window allows the user to adjust the bias calibration 
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process, and to manually input the sensitivity values. Finally, the bottom-right window provides 
a capability to adjust data acquisition parameters and the cut off frequency of the low pass noise 
filters. 

C. Prototype Integration and Test 

Significant efforts addressed integration of the hardware and software components into a 
working end-to-end albeit simplified angular measurement system prototype. At this preliminary 
development stage, a professionally packaged and fully capable system is not the objective, but 
rather to explore feasibility of such a system, to benchmark potential performance, and expose 
any critical problems. Thus, the prototype system is only a temporary arrangement of hardware 
and software devices. The basic integration task consisted of connecting components to one 
another such that communication and data transfer functions between components is reliable and 
accurate. Due to shared laboratory space and equipment, and in particular the angular calibration 
tables and computer workstation, the prototype components required assembly and disassembly 
for each laboratory session. This arrangement was a hindrance to rapid progress. Another 
specific integration problem related to NI-4452 and LabVIEW interfacing. Documentation 
indicated these products are compatible, but after following published and recommended 
guidelines, interfacing was not accomplished. A significant amount of time was devoted to 
resolving this problem. Guidelines were finally abandoned and replaced with the Wyle 
Laboratories core data sampling software heritage. 

After integration tasks were completed, activities focused on prototype testing to 
demonstrate feasibility, assess potential angular measurement performance, and highlight any 
deficiencies. With all components properly connected, one ADXL105 accelerometer and one 
ENC-03J gyro are simultaneously mounted on the angular calibration table. Two methods for 
assessing angular measurement performance are considered. First, the table is rotated manually 
by a specified amount to excite and test the AMS prototype. Comparison of the prototype 
measurement output value with the known table rotation will provide one estimate of the angular 
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accuracy capability. Second, the table is held at rest and measurements are taken by the AMS 
system. Comparisons of the prototype measurement signal stability after bias removal with the 
known zero value for the table will also provide an estimate of the angular measurement 
capability. 

1 0 .0 00 
[deg] 

5 .0 0 0 
0 .0 0 0 
-5 .0 0 0 
-1 0.000 
-1 5 .0 00 

Figure 20 Accelerometer Based AMS Prototype Dynamic Response 
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Figure 21 Accelerometer Based AMS Prototype Static Response 

Measured angular response of the calibration table using the ADXL105 accelerometer 
sensor in the AMS prototype is shown in Figures 20-21. In Figure 20, the table was held at rest 
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initially for approximately 8 s, and then the table was manually rotated to negative, null, positive, 
and back to null orientations with rapid slew motions followed by short holding periods out to 
about 17 s. The angular measurement response in Figure 20 accurately tracks this dynamic 
profile. The negative angular orientation of the table during the rest period was recorded as 12 
deg. The measurement responses appears to be within ±0.1 deg of this known value. Angular 
measurement performance of the AMS prototype using an accelerometer can also be determined 
by observing Figure 21. Figure 21 shows the measurement reading for a period of approximately 
17 s during which the table was held at rest. The response in Figure 21 indicates accelerometer 
signal drift or bias is small and lies approximately within the ±0.1 deg band. Another important 
feature to note from Figure 20 is the high bandwidth of the accelerometer sensor. The ADXL105 
unit shows a rapid response with little overshoot. The small rise time and overshoot features 
visible in Figure 20 most likely represent accurate tracking of physiological limitations during 
manual table excitations. These results solidly confirm the angular measurement accuracy 
projections made in Section II-C based on manufacture’s performance information. 

Measured angular response of the calibration table using the ENC-03J gyro sensor in the 
AMS prototype is shown in Figures 22-23. Recall the ENC-03J gyroscope is mounted on the 
calibration table simultaneously with the ADXL105 accelerometer. Thus, excitations applied to 
the ENC-03J gyro are identical to the excitations applied to the ADXL105 accelerometer. 
Consequently, assuming the accelerometer based measurements are accurate, the angular 
response of ENC-03J based measurement system displayed in Figures 22-23 can be compared 
with the ADXL105 system responses shown in Figures 20-21. Note quantitative comparisons are 
subject to errors in the manually inputed sensitivity values for both sensor units. Figure 22 
indicates the angular measurement response tracks the dynamic motion profile in a rough sense, 
but the gyro output signal drift (indicated by linear growth in angle after integration of the 
constant rate signal) during rest periods is significant. Even worse is the observation that the drift 
is not repeatable and shows strong dependency on response history. Note signal drift becomes 
negative when the gyroscope is turned in the negative direction, and becomes positive when the 
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accuracy capability. Second, the table is held at rest and measurements are taken by the AMS 
system. Comparisons of the prototype measurement signal stability after bias removal with the 
known zero value for the table will also provide an estimate of the angular measurement 
capability. 
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Figure 20 Accelerometer Based AMS Prototype Dynamic Response 
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Figure 21 Accelerometer Based AMS Prototype Static Response 

Measured angular response of the calibration table using the ADXL105 accelerometer 
sensor in the AMS prototype is shown in Figures 20-21. In Figure 20, the table was held at rest 
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All key aspects of the angular measurement system prototype have been investigated. 
Motion detection, data acquisition, data processing, sensor calibration, and attitude construction 
functions were all explored and investigated to a level that confirmed concept feasibility. This 
end-to-end operational capability was tested and demonstrated in the laboratory under real-time 
conditions with actual hardware and software components. Only simplified single-axis 
calibration-construction capability was demonstrated, but this was judged to be sufficient for 
validating feasibility. Angular measurement performance using the accelerometer based 
prototype was confirmed in meeting the design requirements of ±0. 1 deg. Angular measurement 
requirements for the gyro based prototype unfortunately was not confirmed, under the constraints 
of gyro sensor availability. The fundamental problem is poor signal drift characteristics of the 
ENC-03J gyro. Manufacture's performance specifications imply accuracy requirements will be 
difficult to achieve (see Section II-C) unless 1) drift characteristics are highly repeatable and 
stable, 2) precise calibration is achieved, and 3) time efficient test procedures are all employed. 
Methodology for requirement 2 has been offered in Section IV-C and requirement 3 can be 
addressed in future efforts. However, prototype testing confirmed drift characteristics are highly 
erratic and unstable (lack of requirement 1). On the otherhand, if given appropriate MEMS 
performance, the gyro based prototype system should satisfy design objectives. A major 
conclusion of the work is that better performing gyros are required. 

Although the actual prototype system is not professionally packaged, two proposed 
schematics for an AMS prototype specifically tailored for use in the 31 -Inch Hypersonic Tunnel 
are offered. Recall in this application, only one accelerometer package (3 sensors) is needed 
since roll angle is always known. Two options are considered. First, an AMS option using only 3 
accelerometers is considered in Figure 24. Such an AMS has a capability to measure pitch-yaw 
angles to ±0.1 deg. ADXL105 accelerometers are used as representative sensors in Figure 24. 
This AMS arrangement is packaged such that it fits within the tip of typical sting hardware used 
in the 31 -Inch Hypersonic Tunnel. Overall length of the sensor package is about 0.8 in and 
expected mass of the package with mounting plug is approximately 50 g. The second option uses 
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3 accelerometers and 3 gyroscopes and provides a capability to measure pitch-yaw-roll angles to 
±0.1 deg (assuming high performance gyros are used). This option is shown in Figure 25. The 
schematic uses ADXL105 accelerometers and ENC-03J gyros as representative sensors (only for 
dimensional purposes in the gyro case). Expected length of this package is about 1.8 in, and 
expected mass is 100 g with a mounting plug. Section views show that wires can be placed 
within voids beside the accelerometers and gyros. Also, cooling pipes which may be needed in 
the 31 -Inch Hypersonic Tunnel application can be placed in these empty spaces. Utilization of 
motion detection units with multi-axis sensing, such as the MEMS1C units listed in Table 1 or 
other Samsung gyros not considered in this report, can be very useful for reducing weight and 
efficient sensor positioning internal to the sting while maximizing empty space for other 
subsystems. 
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Figure 25 Accelerometer-Gyro AMS Prototype Schematic 
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Section VI 

Conclusions and Recommendations 


Investigation and development of a high accuracy angular measurement capability for 
hypersonic wind tunnel facilities located at NASA LaRC is considered in this effort. 
Specifically, utilization of micro-electro-mechanical sensors including accelerometers and gyros, 
coupled with software driven data acquisition hardware, integrated within a prototype 
measurement system, is addressed. A ±0.1 deg measurement performance in pitch, yaw and roll 
is formulated from wind tunnel constraints and aerodynamic test objectives. A survey of 
commercially available MEMS units is used to project potential measurement performance 
levels that satisfy design requirements. Accelerometer based systems can more easily meet the 
design requirements than gyro based systems. For the gyro based systems to achieve the required 
performance, three factors are needed: 1) highly repeatable and stable drift characteristics, 2) 
highly accurate calibration techniques, and 3) time efficient test procedures. 

A very general analytical framework describing mathematical relationships between multi- 
axis gravitational acceleration and angular rate data to the test article pitch, yaw and roll states is 
developed. An Eulerian attitude description that is consistent with NASA LaRC tunnel facilities 
is selected for use. The framework can address advanced wind tunnel facilities incorporating 
exotic pitch-yaw -roll drive systems, and reduces to simpler applications where roll angle is fixed 
and known (NASA LaRC application). Two accelerometer sensor packages mounted in the sting 
fixed and model fixed reference frames, and one gyro sensor package mounted in the model 
fixed reference frame, are incorporated into the framework. Accelerometer and gyro parameter 
models include mechanical-electrical hardware aspects; specifically bias, sensitivity, and 
azimuth-coning misalignments. Sensor calibration procedures for accelerometers and gyros are 
formulated using nonlinear Gauss-Newton regression logic. These procedures are tested with 
simulated data and without noise and are found to be highly accurate. Angle construction using 
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accelerometers are formulated with both analytic closed-form and iterative Gauss-Newton 
strategies. A two sequence approach using two accelerometer packages facilitates full pitch-yaw- 
roll construction. In the fixed known roll angle case, only one accelerometer package is required. 
Angle construction using gyros is formulated using nonlinear Runge-Kutta numerical 
integration. These procedures are also tested with simulated data without noise and are found to 
be highly accurate. Sensor calibration and attitude construction results validate feasibility of the 
measurement system concept. 

Design of hardware and software components followed by system integration and test 
resulted in a prototype measurement system. Modern software tools are used to create multi- 
functional control instrumentation that reads input signals from a data card, processes and filters 
input signals, calibrates several key sensor parameters, transforms gravitational acceleration and 
angular rate data to angular position states, records data to storage devices, and outputs real-time 
information in numeric and graphical formats. Integration of hardware and software components 
to a laboratory demonstrated end-to-end simplified angular measurement system prototype is 
successfully accomplished. Motion detection, data acquisition, data processing, sensor 
calibration, and attitude construction functions are all explored to a level that confirms concept 
feasibility. Angular measurement performance using the accelerometer based prototype is 
confirmed in meeting the ±0.1 deg design requirement. Angular measurement requirements for 
the gyro based prototype unfortunately is not confirmed. The fundamental problem is poor gyro 
signal drift characteristics. Prototype testing confirmed drift characteristics are highly erratic and 
unstable. However, if given appropriate MEMS performance, the gyro based prototype system 
should satisfy design objectives. A major conclusion of the work is that better performing gyros 
are required. 

The overriding conclusion drawn from this preliminary research effort is that projected 
feasibility for a MEMS based angular measurement system using high grade accelerometers and 
gyros, as well as sophisticated calibration and construction logic, that can satisfy current 
accuracy objectives, is high. However, substantial additional work is necessary to achieve full 
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success. Effects from noise, initial conditions, and sampling rate on the sensor calibration and 
attitude construction algorithms should be addressed and minimized with algorithm 
modification. Rotating earth effects on accelerometer and gyro outputs have been ignored in this 
study and should be addressed in future work. Measurement and control instrumentation 
software needs to be greatly enhanced to achieve a useful working system. General multi-axis 
calibration and construction logic should be incorporated into the software. Alignment 
procedures should also be considered in any future efforts. A broader range of accelerometer and 
gyro MEMS should be requisitioned and sampled for potential performance improvements, 
especially for gyros. Extensive calibration and construction testing including thermal effects 
should also be conducted to build further confidence in the proposed angular measurement 
system. Finally, testing in wind tunnel operations and comparisons with existing measurement 
systems should be pursued. 
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Appendix A 

Statement of Work for PO-L-14635 


Development of a High Accuracy Angular Measurement 
System for LaRC’s Hypersonic Wind Tunnel Facilities 

STATEMENT OF WORK 


1.0 Introduction 

Current angular measurement capability in LaRC’s hypersonic wind tunnel facilities is lacking. 
These current systems are based on primitive mechanical measurement of the model sting pitch 
orientation. Pitch measurement accuracy is insufficient. Yaw and roll orientation measurements 
are not available and must be assumed. Wind tunnel tests are often repeated due to the 
uncertainty of the location of the model after being injected into the tunnel’s test chamber. This 
project focuses on developing a modern and adaptable sensor package capable of measuring 
pitch, roll and yaw while also meeting the size and accuracy requirements stipulated by LaRC’s 
hypersonic wind tunnels. It is believed that a package consisting of gyros, or a combination of 
gyros and accelerometers, will be required. Commercially available MEMS sensor technology is 
projected to have performance levels which can meet both size and accuracy requirements. 

1.1 Scope 

Angular orientation in pitch, roll and yaw can be directly measured by a gyro based sensor 
package. Additionally, pitch and roll orientation can be indirectly measured by an accelerometer 
package with comparison to the gravity vector. Mechanization of angular measurements will be 
based on these methodologies. The contractors shall provide services of a research and 
development nature that follow the specifications below 

2.0 Information provided 

The information related to the angular measurement system will be provided on an as 
needed basis. 

2.1 Materials provided 

LaRC shall provide the following materials 
MEMS sensors 

Laboratory calibration and test equipment 

o Dual axis high precision angle indexing table. The table is built by A.G Davis 
Inc. and is capable of setting 1 -degree increment combinations of pitch and roll 
from 0 to 360 degrees. The table is composed of two separate indexing heads in 
an orthogonal configuration, each accurate to 1 arc second over the full 360- 
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degree angle range. The angular position of each table is computer controlled 
using a serial interface. 

o Single axis high precision indexing table. This table is also computer controlled 
using a serial interface. The table has a range of 360 degrees with a resolution of 
.0001 degrees. 

o Multimeters with scanners for data acquisition such as the HP-3458A 
All hardware required to construct a prototype system 

o Electrical components such as resistors and capacitors 
o Mechanical components such as the housing to hold the sensors 
Computer systems and data acquisition components 
o Computer with data acquisition cards 
Commercial software for developing software for operating the prototype system 
o Visual Basic, Labview, Matlab, ect 
Environment for testing the prototype system 

o Accelerometer laboratory/Wyle Laboratories 
o Temperature chamber 
o Wind tunnel 

Procedures for documenting all work 

3.0 Specifications 

3.1 Collect and analyze wind tunnel angular measurement requirements. 

3.2 Perform a search to find manufacturers of various makes of MEMS gyros and 
accelerometers. Select the best performing sensors based on the manufacturers 
specifications and provide LaRC with a report of the results including all relevant 
information necessary to procure the sensors. 

3.3 Derive the governing equations that model the response of various sensor package 
designs and record the details of each derivation in a document. 

3.4 Evaluate the performance of LaRC procured MEMS gyros and accelerometers by 
conducting tests in LaRC’s accelerometer laboratory and at Wyle Laboratories. 

3.5 Compare test results with tunnel requirements using LaRC approved statistical 
procedures and establish a reasonable goal for both package size and performance. 

3.6 Develop an overall design and operation method for the system. 

3.7 Conduct a preliminary design review and implement suggestions. 

3.8 Assemble a pre-prototype sensor package. 

3.9 Construct a data acquisition system using LaRC supplied hardware . 

3.10 Develop and implement software for calibrating and operating the system. 

3.11 Evaluate the sensor package using LaRC equipment and procedures and compare it’s 
performance with the performance goals. 

3.12 Determine if enhancements or design changes are necessary in order to meet the 
performance goals. If changes are necessary implement the changes and go back to 3.1 1, 
else continue to 3.13. 

3.13 Construct a prototype sensor package. 

3.14 Calibrate and test the prototype system in both LaRC’s accelerometer laboratory and a 
wind tunnel. 

3.15 Collect, analyze and implement feedback based on test results. 
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Appendix B 

Accelerometer Sensor Calibration Listing 


% 


% 3 Axis Accelerometer Package Calibration 

% - Sting Package - 


% % 

% This File Calibrates 3 Accelerometer Oriented in 3 Axis 
% File Name 1 AccCalibOO_Sting 1 

% 

% Research Funded by 

% Measurement & Instrumentation Branch of NASA Langley Research Center 
% NASA LaRC Hampton VA. 23681 
% Instructor in ODU Side: Dr. Newman 
% Program Coded in MATLAB by Si-bok Yu 
% Aerospace Engr. Dept., Old Dominion Univ. 

% Latest Updated in Oct. 2002 


% % 

% 

% Inputs 


% Need Sensor Voltage Output at N Different Pitch, Yaw and Roll Orientations 
% for each of 3 Axis 


% Outputs 

o 

o 

% b % Bias [V] 

% S % Sensitivity [V/g] 

% Om % Coning Angle [rad] 

% A % Azimuth Angle [rad] 

% 

% % 


clear; 

r2d 

d2r 

Tolerance 

G 


= 180/pi; 
= pi/180; 

- le-32; 

= 1 / 


o 

o 


% Load Pitch, Yaw, and Roll Motion, and Voltage Output 
load SimAccData_Sting 


% Initial Guess 

bO 

= 1, 

.23; 

SO 

= 0 , 

,45; 

OmO 

= 0 , 

, 06*d2r; 

AO 

= 0 , 

. 07*d2r; 


Values 

% Bias [V] 

% Sensitivity [V/g] 

% Coning Angle [rad] 

% Azimuth Angle [rad] 


% 
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A_vector = [bO SO OmO AO] ' ; 

nData=length (Phi_m) ; 

for iSen=l:3, 
for ii=l : 100 , 

b = A_vector (1) ; 

S = A_vector (2) ; 

Om = A_vector(3); 

A = A_vector ( 4 ) ; 

for j j®l :nData-l f 

Th = Th_m( j j ) / 

Psi = Psi__m ( j j ) ; 

vx = Vxst ( j j ) ; 
vy = Vyst ( j j ) ; 
vz = Vzst ( j j ) ; 


if iSen == 1 

Z(jj,l)=l; 

Z ( j j , 2 ) = ( “Cos (Om) *cos (Psi) *sin (Th) + sin (Om) *sin (A) *sin ( Psi ) *sin (Th) 

-sin (Om) *cos (A) *cos (Th) 

z ( j j r 3 ) =S* ( sin (Om) *cos (Psi ) *sin (Th) + cos (Om) *sin (A) *sin ( Psi ) *sin (Th) 

-cos (Om) *cos (A) *cos (Th) 

z ( j j , 4) =S* ( sin (Om) *cos (A) *sin (Psi ) *sin (Th) + sin (Om) *sin (A) *cos (Th) 

D ( j j / 1 ) =vx-b 

“S* ( “Cos ( Om ) * cos (Psi) *sin(Th) + sin (Om) *sin (A) *sin ( Psi ) *sin (Th) 

-sin (Om) *cos (A) *cos (Th) 
elseif iSen == 2 
Z(jj,l)=l; 

Z(jj/2)= ( sin (Om) *cos (A) *cos (Psi) *sin (Th) + cos (Om) *sin ( Psi ) *sin (Th) 

+sin (Om) *sin (A) *cos (Th) 

Z ( j j , 3) =S* ( cos (Om) * cos (A) *cos ( Psi ) *sin (Th) - sin (Om) *sin (Psi) *sin (Th) 

+cos (Om) *sin (A) *cos (Th) 

Z(jj,4)=S*( -sin (Om) *sin (A) *cos (Psi) *sin (Th) + sin (Om) *cos (A) *cos (Th) 

D ( j j , 1 ) =vy-b 

-S*( sin (Om) *cos (A) *cos (Psi) *sin (Th) 4- cos (Om) *sin ( Psi ) *sin (Th) 

+sin (Om) *sin (A) *cos (Th) 
elseif iSen — 3 

Z ( j j f 1 ) =1 ; 

z ( j j / 2 ) = ( -sin (Om) *sin (A) *cos (Psi) *sin (Th) - sin (Om) *cos (A) *sin ( Psi ) *sin (Th) 

+cos (Om) *cos (Th) 

Z { j j / 3) =S* ( -cos (Om) *sin(A) *cos (Psi) *sin(Th) - cos (Om) *cos (A) *sin (Psi) *sin (Th) 
-sin (Om) *cos (Th) 

Z ( j j , 4 ) =S* ( -sin (Om) *cos (A) *cos (Psi) *sin (Th) + sin (Om) *sin (A) *sin ( Psi ) *sin (Th) 

D ( j j , 1 ) =vz-b 

-S*( -sin(Om) *sin (A) *cos (Psi ) *sin (Th) - sin (Om) *cos (A) *sin (Psi) *sin (Th) 
+cos (Om) *cos (Th) 
end; 


) *g; 

) *g; 
) *g; 

) *g; 

) *g; 

) *g; 
) *g; 

) *g; 

) *g; 

) *g; 

) *g; 

) *g; 


end; 

dA = inv ( Z' *Z ) * ( Z ' *D); 
A_vector = A_vector + dA; 


error_D = sum(dA. A 2); 
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error_D = sqrt (error_D) ; 
if error_D < Tolerance 
ii 

break; 

end; 

[iSen ii jj/1000 error_D] 
end; 

if iSen == 1 

b_x = A_vector(l); 

S_x = A_vector{2); 

Om_x = A_vector(3); 

A_x = A_vector{4); 
elseif iSen == 2 

b_y = A_vector(l); 

S_y = A_vector(2); 

Om__y - A_vector(3); 

A_y = A_vector ( 4 ) ; 
elseif iSen == 3 

b_z = A_vector(l); 

S_z = A_vector(2); 

Om_z = A_vector (3) ; 

A_z = A_vector(4); 
end; 

end; 

[b_x S_x Om__x*r2d A_x*r2d] 

[b_y S_y Om_y*r2d A_y*r2d] 

[b_z S_z 0m_z*r2d A_z*r2d] 

save AccCalData b_x b_y b_z S_x S_y S_z Om_x 0m_ 


Om_z A_x A_y A_z 
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% % 

% 

% 3 Axis Accelerometer Package Calibration 

% - Model Package - 


% % 

% This File Calibrates .3 Accelerometer Oriented in 3 Axis 
% File Name 1 AccCalibOO_Model ' 

% 

% Research Funded by 

% Measurement & Instrumentation Branch of NASA Langley Research Center 
% NASA LaRC Hampton VA. 23681 
% Instructor in ODU Side: Dr. Newman 
% Program Coded in MATLAB by Si-bok Yu 
% Aerospace Engr. Dept., Old Dominion Univ. 

% Latest Updated in Oct. 2002 
% 

% % 

g 

o 

% Inputs 

% 

% Need Sensor Voltage Output at N Different Pitch, Yaw and Roll Orientations 
% for each of 3 Axis 

g 

o 

g 

o 

% Outputs 

g 

o 

% b % Bias [V] 

% S % Sensitivity [V/g] 

% Om % Coning Angle [rad] 

% A % Azimuth Angle [rad] 

Q, 

O 


clear; 

r2d 

d2r 

Tolerance 

g 


= 180/pi; 
= pi/180; 
= le-32; 

= 1 ; 


% Load Pitch, Yaw, and Roll Motion, and Voltage Output 
load SimAccData Model 


% Initial Guess Values 


bO = 1.23; 

SO = 0.45; 

OmO - 0.-06*d2r; 

AO - 0 . 07*d2r ; 
A_vector = [bO SO OmO 

nData=length ( Phi_m) ; 

for iSen=l:3, 
for ii=l:100, 

b = A_vector(l 


% Bias [V] 

% Sensitivity [V/g] 

% Coning Angle [rad] 

% Azimuth Angle [rad] 
AO] ' ; 


% 


% 
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S = A_vector (2 ) ; 
Om = A_vector(3); 
A = A__vector (4) ; 

for j j=l :nData-l, 


Phi = Phi m ( j j ) ; 

Th 

= Th_m ( j j ) ; 

Psi 

vx = Vxm ( j j ) ; 

vy 

= Vym ( j j ) ; 

vz 


Psi_m( j j ) ; 

Vzm ( j j ) ; 

if iSen == 1 
Z(jj,l)=l; 

Z(jj,2)= ( -cos (Om) *cos (Psi) *sin (Th) 

+ (sin (Om) *sin (A) *cos (Phi) +sin (Om) *cos (A) *sin (Phi) ) *sin (Psi) *sin (Th) 

+ (sin (Om) *sin (A) *sin (Phi) -sin (Om) *cos (A) *cos (Phi) ) *cos (Th) ) 

Z ( j j , 3) =S* ( sin (Om) *cos (Psi) *sin (Th) 

+ (cos (Om) *sin (A) *cos (Phi) +cos (Om) *cos (A) *sin (Phi) ) *sin (Psi) *sin (Th) 

+ (cos (Om) *sin (A) *sin (Phi) -cos (Om) *cos (A) *cos (Phi) ) *cos (Th) 

Z ( j j, 4)=S* ( sin (Om) *cos (A) *cos (Phi ) -sin (Om) *sin (A) *sin (Phi) ) *sin(Psi) *sin(Th) 

+ (sin (Om) *sin (A) *cos (Phi) +sin (Om) *cos (A) *sin (Phi) ) *cos (Th) 

D ( j j , 1 ) =vx-b 

-S* ( -cos (Om) *cos (Psi) *sin (Th) 

+ (sin (Om) *sin (A) *cos (Phi) +sin (Om) *cos (A) *sin (Phi) ) *sin (Psi) *sin (Th) 

+ (sin (Om) *sin (A) *sin (Phi) -sin (Om) *cos (A) *cos (Phi) ) *cos (Th) 
elseif iSen == 2 
Z ( j j » 1) =1; 

Z ( j j » 2) = { sin (Om) *cos (A) *cos (Psi) *sin (Th) 

+ (cos (Om) *cos (Phi) sin (Om) *sin (A) *sin (Phi) ) *sin (Psi) *sin (Th) 

+ (cos (Om) *sin (Phi) +sin (Om) *sin (A) *cos (Phi) ) *cos (Th) 

Z ( j j , 3)=S* ( cos (Om) *cos (A) *cos (Psi) *sin (Th) 

+ (-sin (Om) *cos (Phi) -cos (Om) *sin (A) *sin (Phi) ) *sin (Psi) *sin (Th) 

+ (-sin (Om) *sin (Phi) +cos (Om) *sin (A) *cos (Phi) ) *cos (Th) 

Z (j j,4)=S* ( -sin (Om) *sin(A) *cos (Psi) *sin(Th) 

-sin (Om) *cos (A) * sin (Phi) *sin (Psi) *sin (Th) 

+sin (Om) *cos (A) *cos (Phi) *cos (Th) 

D ( j j , 1) =vy-b 

-S* ( sin (Om) *cos (A) *cos (Psi) *sin (Th) 

+ (cos (Om) *cos (Phi) - sin (Om) *sin (A) *sin (Phi) ) *sin (Psi) *sin (Th) 

+ (cos (Om) *sin (Phi) +sin (Om) *sin (A) *cos (Phi) ) *cos (Th) 
elseif iSen == 3 
Z ( j j / 1) =1; 

Z ( j j / 2 ) = ( -sin (Om) *sin (A) *cos (Psi) *sin (Th) 

+ (-sin (Om) *cos (A) *cos (Phi) -cos (Om) * sin (Phi) ) *sin ( Psi ) *sin (Th) 

+ (-sin (Om) *cos (A) *sin (Phi) +cos (Om) *cos (Phi) ) +cos (Th) 

Z (j j, 3) =S* ( -cos (Om) *sin (A) *cos (Psi) *sin (Th) 

+ (-cos (Om) * cos (A) *cos (Phi) +sin (Om) *sin (Phi) ) *sin (Psi) *sin (Th) 

+ (-cos (Om) *cos (A) *sin (Phi) -sin (Om) *cos (Phi) ) *cos (Th) 


: g; 


Z(jj, 

4) 

=s* 

{ -sin(Om) 

*cos (A) 

*cos (Psi) 

*sin (Th) 




+sin (Om) 

+ sin (A) 

*cos (Phi) 

*sin ( Psi ) 




+sin (Om) 

*sin (A) 

*sin (Phi) 

*cos (Th) 

D ( j j , 

1) 

=vz- 

-b 






-S* ' 

( -sin(Om) 

*sin (A) 

*cos (Psi) 

*sin (Th) 


k sin (Th) 


+ (-sin (Om) *cos (A) *cos (Phi) -cos (Om) *sin (Phi) ) *sin(Psi) *sin (Th) 
+ (-sin (Om) *cos (A) *sin (Phi) +cos (Om) *cos (Phi) ) *cos (Th) 
end; 


) *g; 
) *g; 


) *g; 


) *g; 


k g; 


) *g; 


) *g; 


) *g; 


) *g; 


k g; 


) *g; 


end ; 
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dA = inv { Z 1 *Z ) * ( Z 1 *D ) ; 

A_ vec t° r = A_vector + dA; 

error_D = sum(dA. A 2); 
error_D = sqrt (error_D) ; 
if error_D < Tolerance 
ii 

break; 

end; 

[iSen ii jj/1000 error_D] 
end; 

if iSen == 1 

b_x2 = A_vector(l); 

S_x2 = A_vector(2); 

0m_x2 = A_vector(3); 

A_x2 = A_vector(4); 
elseif iSen == 2 

b_y2 = A_vector(l); 

S__y2 = A_vector(2); 

0m_y2 = A_vector(3); 

A_y2 = A_vector(4); 
elseif iSen == 3 

b_z2 = A_vector(l); 

S_z2 = A_vector(2); 

0m_z2 = A_vector(3); 

A_z2 = A_vector(4); 
end; 

end; 

[b_x2 S_x2 Om_x2*r2d A_x2*r2d] 

[b_y2 S_y2 Om_y2*r2d A_y2*r2d] 

[b_z2 S_z2 Om_z2*r2d A_z2*r2d] 

save AccCalData2 b_x2 b_y2 b_z2 S_x2 S_y2 S_z2 0m_x2 0m_y2 0m_z2 A_x2 A_y2 
A z2 


I 

i 



91 


Appendix C 

Accelerometer Attitude Construction Listing 


% 


% Pitch, Yaw, and Roll Angle Calculation from Accelerometer Package 

% - Sting & Model Packages - 

% 

% % 

% This File Calibrates Pitch, Yaw, and Roll Angle from Acc. Oriented in 3 Axis 
% File Name , Acc_Angle* 

% 

% Research Funded by 

% Measurement & Instrumentation Branch of NASA Langley Research Center 
% NASA LaRC Hampton VA. 23681 
% Instructor in ODU Side: Dr. Newman 
% Program Coded in MATLAB by Si-bok Yu 
% Aerospace Engr. Dept., Old Dominion Univ. 

% Oct. 2002 


% 

% This File Takes Accelerometer Package Calibration Information from 
% 1 AccCalibO Ousting ' & 1 AccCalibOO_Model 1 
% The Information Includes 
% b % Bias [V] 

% S % Sensitivity [V/g] 

% Om % Coning Angle [rad] 

% A % Azimuth Angle [rad] 

% 

% Outputs 

% 

% This File Gives Output 
% Pitch Angle [rad] 

% Yaw Angle [rad] 

% These Outputs are used as Inputs to Calculate the Additional Output 
% Roll Angle [rad] 

o 

"6 

% 


clear; 

r2d 

d2r 

Tolerance 

g 


= 180/pi; 
= pi/180; 
= le-10 ; 

= l; 


o 

"O 

% Load Calibration Information for Sting Package 
% Load Bias, Sensitivity, Coning Angle, and Azimuth Angle Data 
load AccCalData 

% Load Calibration Information for Model Package 
load AccCalData2 


a 

o 
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% Load Voltage Output from the Sting Package 
% at Corresponding Pitch, Yaw, and Roll Motion 
load SimAccData2_Sting 

% Load Voltage Output from the Model Package 
load SimAccData2_Model 

% Initial Value for Computing Orientation of the Model 


PhiO 

= 0.1; 

ThO 

= 0.1; 

PsiO 

= 0.1; 

nData 

= length (Vxst) 


% 

% 

% Find Angles Using Iterative Method (Gauss Newton) 

% 

% 

for jj=l:nData, 

Vx = Vxst ( j j ) ; 

Vy = Vyst ( j j ) ; 

Vz = Vzst ( j j ) ; 

A_vector = [ThO PsiO] f ; 

for ii=l:10, 

Th = A_vector(l); 

Psi = A_vector(2); 

Z (1, 1) =S_x* ( -cos (Om_x) *cos (Psi) *cos (Th) 

+sin (Om_x) *sin (A_x) * sin (Psi) *cos (Th) 
+sin (Om_x) *cos (A_x) *sin (Th) 

Z ( 1, 2 ) =S__x* ( +cos (Om_x) * sin (Psi ) *sin (Th) 

+sin (Om_x) *sin (A_x) *cos (Psi) *sin (Th) 

Z (2, 1) =S_y* ( + sin (Om_y) *cos (A_y) *cos (Psi) *cos (Th) 

+ cos (Om_y) * sin (Psi ) *cos (Th) 

-sin (Om_y) *sin (A_y) *sin (Th) 

Z (2, 2 ) =S_y* ( -sin (Om_y) *cos (A y ) *sin ( Psi ) *sin (Th) 

+cos (Om_y) *cos (Psi ) *sin (Th) 

Z (3, 1) =S_z* ( -sin (Om_z) *sin (A_z) *cos (Psi) *cos (Th) 
-sin (Om_z) *cos (A_z) *sin (Psi) *cos (Th) 
-cos (Om_z) *sin (Th) 

Z (3, 2) =S_z* ( +sin (Om__z) *sin (A_z) * sin (Psi) *sin (Th) 
-sin(Om z)*cos(A z ) *cos ( Psi ) *sin (Th) 


D ( 1 ) =Vx-b_x 

i 

CO 

1 

X 

* 

-cos (Om x) *cos ( Psi ) *sin (Th) 

+sin (Om_x) *sin (A_x) *sin(Psi)*sin(Th) 
-sin (Om_x) *cos (A_x) *cos (Th) 

) *g; 

D ( 2 ) =Vy-b_y 

-S_y* ( 

sin (Om_y) *cos (A_y) *cos (Psi) *sin (Th) 
+cos (Om_y) *sin ( Psi ) *sin (Th) 

+sin (Om y)*sin(A y) *cos (Th) 

) *g; 

D ( 3) =Vz-b_z 

-S_z * ( 

-sin (Om_z) *sin (A__z) *cos (Psi) *sin (Th) 
-sin (Om_z) *cos (A_z) * sin (Psi) *sin (Th) 
+cos (Om z)*cos(Th) 

) *g; 


) *g; 
) *g; 

) *g; 
) *g; 

) *g; 
) *g; 


Q. 

"O 


93 


dA = inv ( Z* *Z ) * ( Z 1 *D' ) ; 
A_vector = A_vector + dA; 

error_D = sum(dA. A 2); 
error_D - sqrt (error_D) ; 
if error_D < Tolerance 
[jj ii error_D] 

break; 

end; 

end; % ii Loop 

Vx2 = Vxm C j j ) ; 

Vy2 = Vym ( j j ) ; 

Vz2 = Vzm ( j j ) ; 

Th = A_vector ( 1 ) ; 

Psi = A__vector (2 ) ; 

A_vector2 = [ Phi 0 ] 1 ; 

for ii=l:10. 

Phi = A vector2(l); 


Z(l,l)=S_x2*( 

+ (-sin (0m_x2 ) *sin (A_x2 ) *sin ( Phi ) +sin (Om_x2 ) *cos (A__x2 ) *cos ( Phi) )*sin(Psi)*sin( 
+ (sin (0m_x2) *sin (A_x2) *cos (Phi) +sin (0m_x2) *cos (A_x2) *sin (Phi) ) *cos (Th) ) *g; 

Z (2 , 1 ) =S_y2* ( 

+ (-cos (0m_y2) *sin (Phi) -sin (0m__y2) *sin (A_y2) *cos (Phi) ) * sin (Psi) *sin (Th) 

+ (cos (0m_y2) *cos (Phi) -sin (0m_y2) *sin (A_y2) *sin (Phi) ) *cos (Th) ) *g; 

Z ( 3, 1 ) =S_z2* ( 

+ (sin (0m_z2) *cos (A_z2) *sin (Phi) -cos (0m_z2) *cos (Phi) ) *sin (Psi) *sin (Th) 

+ (-sin (0m_z2 ) *cos (A_z2 ) *cos ( Phi) -cos (0m_z2 ) *sin (Phi ) ) *cos (Th) ) *g; 

D (1) =Vx2-b_x2 -S_x2* ( -cos (Om_x2) *cos (Psi) *sin (Th) 

+ (sin (Om_x2 ) *sin (A_x2 ) *cos ( Phi ) +sin (0m_x2 ) *cos (A_x2 ) *sin ( Phi) ) *sin(Psi) *sin (Th) 
+ ( sin (0in_x2 ) *sin (A__x2 ) *sin ( Phi ) - sin (Om_x2 ) *cos (A__x2 ) *cos ( Phi ) ) *cos (Th) ) *g; 

D (2) =Vy2-b_y2 -S_y2* ( sin (0m_y2) *cos (A_y2 ) *cos ( Psi ) *sin (Th) 

+ ( cos (Om_y2 ) *cos { Phi ) -sin (Om_y2 ) *sin (A_y2 ) *sin (Phi ) ) * sin (Psi ) *sin (Th) 

+ (cos (Om_y2) *sin (Phi) +sin (Om_y2) *sin (A_y2) *cos (Phi) ) *cos (Th) ) *g; 

D ( 3 ) =Vz2-b_z2 -S_z2* ( -sin(0m_z2) *sin ( A_z2 ) *cos ( Psi ) *sin (Th) 

+ (-sin (0m_z2 ) *cos (A_z2 ) *cos (Phi) -cos (0m_z2 ) *sin ( Phi ) ) *sin ( Psi ) *sin (Th) 

+ (-sin (0m_z2 ) *cos (A__z2 ) *sin ( Phi) +cos (Om_z2 ) *cos ( Phi ) ) *cos (Th) ) *g; 


dA2 = inv ( Z 1 *Z ) * ( Z 1 *D 1 ) ; 
A_vector2 - A_vector2 + dA2 ; 

error^D = sum (dA2 . ^2 ) ; 
error^D = sqrt (error_D) ; 
if error_D < Tolerance 
[jj ii error_D] 

break; 

end; 

end; % ii Loop 
Phi = A vector2(l); 
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Phi_n(jj) = Phi; 
Th_n(jj) = Th; 
Psi_n(jj) = Psi; 

end; % jj loop 


save AccAngle 


o _ ... _ _ _ „ _ _ _ ____________ 

0 — — — . — 

o 

o 

% Closed-Form Solution 

o 

o 

o ____ ________ 

"O — 


o 


% Bias and Sensitivity Matrices for Sting Package 
B_sn_st = [b_x b_y b_z] . 1 ; 

S_sn_st = [ S_x 0 0 

0 S_y 0 

0 0 S_z ] ; 

% Bias and Sensitivity Matrices for Model Package 
B_sn_m = [b__x2 b_y2 b_z2] . f ; 

S_sn_m = [S_x2 0 0 

0 S_y2 0 

0 0 S_z 2 ] ; 


for qq = l:nData 


Vx_st = Vxst(qq); 
Vy_st = Vyst(qq); 
Vz__st = Vzst(qq); 


Vx_m = Vxm(qq); 
Vy_m = Vym(qq); 
Vz_m = Vzm(qq); 


V_sn_st = [ Vx_st Vy__st Vz_st]. f ; 
V_sn_m = [ Vx_m Vy_m Vz_m ] . ’ ; 


T_sn_x_st_st 1 
T_sn_y_st_st2 
T sn z st st3 


[ cos (Om_x) 

[-sin (Om_y) *cos (A_y) 
[ sin (Om_z ) *sin (A_z ) 


sin (Om_x) *sin (A__x) 
cos (Om_y) 

-sin(Om z)*cos(A__z) 


-sin (Om_x) *cos (A_x) ] ; 
sin (Om_y) *sin (A_y) ] ; 
cos (Om z ) ] ; 


T_sn_st_st = [T_sn_x_st_st 1 
T_sn_y_st_st2 
T_sn_z_st_st 3 ] ; 


T_sn_x_m_ml 
T_sn_y_m_m2 
T sn z m m3 


[ cos(0m_x2) 

[-sin (0m_y2) *cos (A_y2) 
[ sin (0m_z2 ) *sin (A z2 ) 


sin (0m_x2 ) *sin (A_x2 ) 
cos (0m_y2 ) 

-sin(Om z2)*cos(A z2) 


-sin (0m_x2) *cos (A_x2) 
sin (0m_y2) *sin (A_y2) 
cos (Om z2 ) ] ; 
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T_sn_m_m = [T_sn__x_m_ml 
T__sn__y_m__m2 
T_sn__z_m_jn3 ] ; 

G_st = inv (T_sn_st_st) *inv (S_sn_st ) * (V_sn_st - B_sn_st) ; 

G_m = inv (T_sn_m_m) *inv (S_sn_m) * (V_sn_m - B_sn_m) ; 

G_x_st = G_st (1) / 

G y st — G_st ( 2 ) / 

G~z~st = G~st(3); 

G_x_m = G_m ( 1 ) / 

G_y_m = G_m (2) ; 

G__z_m = G_m (3) ; 

psi{qq) = atan2 (G_y_st , -G_x__st ) ; 

thta (qq) = at an (sqrt (G__x_st A 2+G_y_st ^2 ) /G_z_st) *sign (G_x_st /1/cos (psi (qq) ) ) ; 
phi(qq) = atan2 (G__y_m*G_z_st-G_z_m*G_y_st , G_y_m*G_y_st +G_z_m*G_z_st ) ; 

end % End of pp Loop 
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Appendix D 

Gyro Sensor Calibration Listing 


3 Axis Gyroscope Package Calibration 


% % 

% This File Calibrates 3 Gyroscope Oriented in 3 Axis 
% File Name ' GyroCalibOO 1 

% 

% Research Funded by 

% Measurement & Instrumentation Branch of NASA Langley Research Center 
% NASA LaRC Hampton VA. 23681 
% Instructor in ODU Side: Dr. Newman 
% Program Coded in MATLAB by Si-bok Yu 
% Aerospace Engr. Dept., Old Dominion Univ. 

% Oct. 2002 


% % 

% 

% Inputs 

% 


% Need Sensor Voltage Output at N Different Yaw/Roll Orientations 
% and Pitch/Yaw/Roll Rates for each of 3 Axis 

% 

% 

% Outputs 

% 

% b % Bias [V] 

% S % Sensitivity [V/g] 

% Om % Coning Angle [rad] 

% A % Azimuth Angle [rad] 

% 

% % 


clear; 
r2d= 180/pi; 
d2r=pi/180; 
Tolerance=le-10; 


% Load Pitch, Yaw, and Roll Motion, and Voltage Output 
% at the Corresponding Pitch, Yaw, and Roll Motion 
load SimData 


% % 

% Initial Guess Values 
bO = 0.05; % Bias [V] 

SO = 1.25; % Sensitivity [V/g] 

OmO = 0.1*d2r; % Coning Angle [rad] 

AO = 0.1*d2r; % Azimuth Angle [rad] 

A vector = [bO SO OmO AO] '; 
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nData^length ( Phi_m) ; 

for iSen=l:3, 
for ii=l:20, 

b = A_vector (1) ; 

S = A_vector (2) ; 

Om = A_vector(3); 

A = A__vector ( 4 ) ; 

for j j =1 : nData-1, 

Phi = Phi__m ( j j ) ; 

Th = Thjn(jj) ; 

Psi = Psi_m ( j j ) ; 

dPhi = dPhi_m(jj); 
dTh = dTh_m(jj); 
dPsi = dPsi_m( j j ) ; 


vx=Vx_m ( j j ) ; 
vy=Vy_m(j j) ; 
vz=Vz_m { j j ) ; 


if iSen == 1 


z ( j j , 1 ) = 1 ; 

Z(jj,2)= 

cos (Om) * (sin (Psi) *dTh+dPhi) 

+sin (Om) * sin (A) * (cos (Phi) *cos (Psi) *dTh+sin (Phi) *dPsi) 
-sin (Om) *cos (A) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi ) *dPsi) ; 


Z(jj,3)=S* 

( -sin (Om) * (sin (Psi) *dTh+dPhi) 

+cos (Om) *sin (A) * (cos (Phi) *cos (Psi) *dTh+sin (Phi) *dPsi) 
-cos (Om) * cos (A) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi) *dPsi) 

) ; 

Z ( j j , 4 ) =S* 

( sin (Om) *cos (A) * (cos (Phi) *cos (Psi) *dTh+sin (Phi) *dPsi) 
+sin (Om) *sin (A) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi) *dPsi) 

) ; 

D(j j, 1 ) =vx- 

-b 


-S* 

( cos (Om) * (sin (Psi) *dTh+dPhi) 

+sin (Om) *sin (A) * (cos (Phi) *cos (Psi) *dTh+sin (Phi) *dPsi) 
-sin (Om) *cos (A) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi) *dPsi) 

) ; 

elseif iSen == 2 


z ( j j , 1 ) =1 ; 
Z(jj,2) = 

-sin (Om) *cos (A) * (sin (Psi) *dTh+dPhi) 

+cos (Om) * (cos (Phi) *cos (Psi) *dTh+sin (Phi) *dPsi) 

+sin (Om) *sin (A) * (-sin(Phi) *cos (Psi) *dTh+cos (Phi) *dPsi) ; 


Z ( j j / 3) =S* i 

( -cos (Om) *cos (A) * (sin (Psi ) *dTh+dPhi) 

-sin (Om) * (cos (Phi) *cos (Psi) *dTh+sin (Phi) *dPsi) 

+cos (Om) *sin (A) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi ) *dPsi) 

) ; 

Z ( j j , 4 ) =S* - 

( sin (Om) *sin(A)*(sin(Psi) *dTh+dPhi) 

+ sin (Om) *cos (A) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi) *dPsi) 

) ; 

D(jj, l)=vy- 

-b 


-S* 1 

; -sin (Om) *cos (A) * (sin (Psi) *dTh+dPhi) 

+cos (Om) * (cos (Phi) *cos ( Psi ) *dTh+sin (Phi) *dPsi) 

+sin (Om) *sin (A) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi) *dPsi) 

) ; 


elseif iSen == 3 

z ( j j / 1 ) = l ; 

Z ( j j , 2) = sin (Om) *sin (A) * (sin(Psi) *dTh+dPhi) 


-sin (Om) + cos (A) * (cos (Phi) *cos (Psi) *dTh+sin ( Phi ) *dPsi) 
+COS (Om) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi) *dPsi) ; 
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) 


Z ( j j / 3) =S* ( cos (Om) *sin (A) * (sin (Psi) *dTh+dPhi) 

-cos (Om) *cos (A) * (cos (Phi) *cos (Psi) *dTh+sin (Phi) *dPsi) 

-sin (Om) * (-sin (Phi) * cos (Psi) *dTh+cos (Phi) *dPsi) 

Z ( j j / 4 ) =S* ( sin (Om) *cos (A) * (sin ( Psi ) *dTh+dPhi ) 

+sin (Om) *sin (A) * (cos (Phi) *cos (Psi) *dTh+sin (Phi) *dPsi) ) 

D ( j j , 1 ) =vz-b 

— S * ( sin(Om) *sin (A) * (sin(Psi) *dTh+dPhi) 

-sin (Om) *cos (A) * (cos ( Phi ) *cos ( Psi ) *dTh+sin ( Phi ) *dPsi ) 

+cos (Om) * (-sin (Phi) *cos (Psi) *dTh+cos (Phi) *dPsi) ) 

end; 

end; 

dA = inv ( Z 1 * Z ) * ( Z 1 * D ) ; 

A_vector = A_vector + dA; 

%if A_vector(4) > 2*pi 
% A_vector(4) = A_vector(4) - 2*pi; 

%end; 

error_D = sum(dA. A 2); 
error_D = sqrt (error_D) ; 
if error_D < Tolerance 
[iSen ii] 
break; 
end; 

end; 

if iSen == 1 

b_x = A_vector(l); 

S_x = A_vector(2); 

Om_x = A_vector(3); 

A_x = A_vector ( 4 ) ; 
elseif iSen == 2 

b_y = A_vector(l); 

S_y = A_vector{2); 

Om_y = A_vector(3); 

A_y = A_vector(4); 
elseif iSen == 3 

b_z = A_vector(l); 

S_z = A_vector(2); 

Om_z = A_vector(3); 

A_z = A_vector ( 4 ) ; 
end; 

end; 

[b_x S_x Om_x*r2d A_x*r2d] 

[b_y S_y Om_y*r2d A_y*r2d] 

[b_z S_z Om_z*r2d A_z*r2d] 

save GyroCalData b_x b_y b_z S_x S_y S_z Om_x Om_y Om_z A_x A_y A_z 
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Appendix E 

Gyro Attitude Construction Listing 


% % 

% 

% Pitch, Yaw, and Roll Angle Calculation from Gyroscope Package 

% 

% % 

% This File Calibrates Pitch, Yaw, and Roll Angle from Gyroscopes Oriented in 3 Axi 
% File Name 1 Gyro_Angle' 

% 

% Research Funded by 

% Measurement & Instrumentation Branch of NASA Langley Research Center 
% NASA LaRC Hampton VA. 23681 
% Instructor in ODU Side: Dr. Newman 
% Program Coded in MATLAB by Si-bok Yu 
% Aerospace Engr. Dept., Old Dominion Univ. 

% Feb. 2002 
% 

% % 

% 

% This File Takes Gyroscope Package Calibration Information from 
% ’GyroCalibOO ? 

% The Information Includes 
% b % Bias [V] 

% S % Sensitivity [V/g] 

% Om % Coning Angle [rad] 

% A % Azimuth Angle [rad] 

% 

% Outputs 

% 

% This File Gives Output 
% Pitch Rate [rad/s] 

% Yaw Rate [ra d/s] 

% Roll Rate [rad/s] 

% which are Integrated To Yield 
% Pitch Angle [rad] 

% Yaw Angle [rad] 

% Roll Angle [rad] 

% 

% % 


clear; 
r2d=180/pi ; 
d2r=pi/180 ; 
Tolerance=le-10 ; 


% Load Calibration Information 

% Load Bias, Sensitivity, Coning Angle, and Azimuth Angle Data 
load GyroCalData 

% % 
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% Load Voltage Output from the Gyroscope Package 
% at the Corresponding Pitch, Yaw, and Roll Motion 
load SimData2 


% Initial Orientation of the Model 
Phi0=0 ; 

ThO =0; 

Psi0=0 ; 

% 

% Initial Guess Values 

P0= (Vx_m ( 1 ) - b_x) /S_x; Q0=(Vy_m(l) - b_y) /S_y; R0=(Vz_m(l) - b_z) /S_z ; 


Tmt -[1 

sin ( PsiO ) 

0; 

0 

cos (PhiO) *cos (PsiO) 

sin (PhiO) ; 

0 

-sin (Phi 0) *cos (PsiO) 

cos (PhiO) ] ; 

dEulerO 

= inv(Tmt) * [P0 Q0 R0] 

t . 

r 


dPhiO=dEulerO (1) ; dThO =dEulerO (2) ; dPsiO^dEulerO (3) ; 
nData=length (Vx__m) ; 

X= [ 0 0 0 ] 1 ; 


% Find Angles with Iteration Method 
for jj=l:nData, 

V_x=Vx_m ( j j ) ; 

V_y=Vy_m (jj); 

V_ z= Vz_m ( j j ) ; 

save GyroEOM_Inputs V_x V_y V_z 


% 


Use Runge-Kutta 5th Order Method from ’Numerical Method for Engineering, ' 


Chapra & Canale, McGraw Hill, P. 604 
kl=GyroEOM_rk5 ( X ( : , j j ) 
k2=GyroEOM_rk5 { X ( : , j j ) +l/4*dt*kl ! 
k3=GyroEOM_rk5 ( X(:,jj) +l/8*dt*kl ! 
k4=GyroEOM_rk5 ( X ( : , j j ) -l/2*dt*k2 1 
k5=GyroEOM_rk5 ( X ( : , j j ) +3/16*dt*kl 1 
k6=GyroEOM_rk5 ( X ( : , j j ) -3/7*dt*kl' 

+ 8/7 *dt *k5 1 


1990 

/ 1 ( j j ) ) ; 

, t ( j j ) +l/4*dt) ; 

+l/8*dt*k2' ,t(jj)+l/4*dt); 

+dt *k3 1 , t (j j)+l/2*dt) ; 

+9/16*dt*k4 ? , t ( j j ) +3/4 *dt ) ; 

+2/7*dt*k2 ' +12/7*dt*k3 1 -12/7*dt*k4' 

, t ( j j ) +dt ) ; 


X ( : f j j tl ) =X ( : , j j ) + (1/90 * (7*kl + 32*k3 + 12*k4 + 32*k5+7*k6) ) 1 *dt; 
clear kl k2 k3 k4 k5 k6; 


[jj X ( 1 , j j ) *r2d X (2, j j ) *r2d X(3,jj)*r2d] 

Phi_n ( j j ) = X ( 1 , j j > ; Th_n(jj) = X ( 2 , j j ) ; Psi_n(jj) == X (3, j j ) ; 
end; % jj loop 
save GyroAngle 
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o\° o\° M o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° 


function f=GyroE0M_rk5 (X, t) 

a _ 

"6 

% 

% Equation of Motion of the Gyroscope Sensor Package 

% 

% 

% This is in the Form of Non-Linear State-Space Equation 
% This File is being Called by the Main File 'Gyro_Angle' 

% 

% Research Funded by 

% Measurement & Instrumentation Branch of NASA Langley Research Center 
% NASA LaRC Hampton VA. 23681 
% Instructor in ODU Side: Dr. Newman 
% Program Coded in MATLAB by Si-bok Yu 
% Aerospace Engr. Dept., Old Dominion Univ. 

% Mar. 2002 

o 

■?> 

% 

% 

% This File Takes Accelerometer Package Calibration Information from 
% ’ GyroCalibOO ' 

% The Information Includes 
% b % Bias [V] 

% S % Sensitivity [V/g] 

% Om % Coning Angle [rad] 

A % Azimuth Angle [rad] 

Inputs 

( Gyroscope Package Voltage Output ) 

Vx, Vy, Vz 


Load Calibration Information from 1 GyroCalibOO 1 
Load Bias, Sensitivity, Coning Angle, and Azimuth Angle Data 
oad GyroCalData 


Load Voltage Calibration Information from ’GyroCalibOO 1 
load GyroEOM_Inputs 

Phi =X(1); 

Th =X ( 2 ) ; 

Psi =X ( 3 ) ; 


Q. 

o 


% 


% 


% 


% 


% 


f (1) = ( (cos (Om_y) *cos (Om_z) +sin (Om_y) *sin (A y ) *sin (Om_z) *cos (A_z) ) 

/ (cos ( Om_x ) *cos (Om_y) *cos (Om_z) +cos (Om_x) *sin (Om_y) *sin (A_y) 

*sin (Om_z) *cos (A_z) + sin (Om_y) *cos (A_y) *sin (Om_x) *sin (A_x) *cos (Om_z) 
-sin (Om_y) *cos (A_y ) *sin (Om_x) *cos (A_x) *sin (Om_z) *cos (A_z ) + sin (Om_z ) 
*sin (A_z) *sin (Om_x) *sin (A__x) *sin (Om_y) *sin (A_y) +sin (Om_z) *sin (A_z) 
*sin (Om_x) *cos (A_x) *cos (Om_y ) ) -sin ( Psi ) *cos ( Phi ) / cos ( Psi ) / (cos ( Phi ) ^2 
+ sin (Phi) ^2) *sin (Om__y) * (cos (A_y) *cos (Om_z) +sin (A_y) *sin (Om_z) 

*sin (A_z ) ) / (cos (Om_x) *cos (Om_y) *cos (Om_z ) +cos (Om_x) *sin (Om_y) 

*sin (A_y) *sin (Om_z ) *cos ( A_z ) +sin (Om_y) *cos (A_y) *sin (Om_x) *sin (A_x) 
*cos (Om_z ) -sin (Om_y) *cos (A__y ) *sin (Om_x) *cos (A_x) *sin (Om_z ) *cos (A_z ) 
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+sin (Om_z) *sin (A_z) *sin (Om_x) ‘sin (A_x) *sin (Om_y) ‘sin (A_y) +sin (Om_z) 
*sin (A_z ) ‘sin (Om_x) ‘cos (A_x) ‘cos (Om_y) ) -sin (Psi) ‘sin (Phi) /cos (Psi) 

/ (cos (Phi) A 2+sin (Phi) A 2) *sin (Om_z) * (-sin (Om_y) ‘cos (A_y) *cos (A z) 

+cos ( Om_y ) *sin (A_z) ) / (cos (Om_x) *cos (Om_y) ‘cos (Om_z) +cos (Om_x) 

*sin (Om_y) *sin (A_y) ‘sin (Om_z) *cos (A_z) +sin (Om_y) ‘cos (A_y) ‘sin (Om_x) 
*sin (A_x) * cos (Om_z) -sin (Om_y) ‘cos (A_y) ‘sin (Om_x) *cos (A_x) ‘sin (OitTz) 
*cos (A_z) +sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin (Om_y) ‘sin (A_y) 
+sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) ) /S_x* (V_x-b_x) 

+ (-sin (Om_x) * (sin (A_x) ‘cos (Om_z) -cos (A_x) *sin (Om_z) *cos (A_z) ) 

/ (cos (Om_x) *cos (Om_y) ‘cos (Om_z) +cos (Om_x) *sin (Om_y) ‘sin (A_y) 

*sin (Om_z) *cos (A_z) +sin (Om_y) ‘cos (A_y) ‘sin (Om_x) *sin (A_x) ‘cos (Om_z) 
-sin (Om_y) *cos (A_y) ‘sin (Om_x) ‘cos (A _x) ‘sin (Om_z) *cos (A_z) +sin (Om_z) 
*sin (A_z) ‘sin (Om_x) *sin (A_x) ‘sin (Om_y) *sin (A_y) +sin (Om_z) ‘sin (A_z) 
*sin (Om_x) *cos (A_x) *cos (Om_y) ) -sin ( Psi ) *cos (Phi) /cos (Psi) / (cos (Phi) A 2 
+sin (Phi ) A 2) * (cos (Om_x) ‘cos (Om_z) +sin (Om_x) *cos (A_x) ‘sin (Om_z) 

*sin (A_z) ) / (cos (Om_x) *cos (Om_y) *cos (Om_z) +cos (Om_x) ‘sin (Om_y) 

*sin (A_y) *sin(Om_z) *cos (A_z) +sin (Om_y) ‘cos (A_y) *sin (Om_x) ‘sin (A_x) 
‘cos (Om_z) -sin (Om_y) *cos (A_y) ‘sin (Om_x) *cos (A _x) ‘sin (Om_z) *cos (A_z) 
+sin (Om_z) ‘sin (A_z) ‘sin (Om_x) ‘sin (A_x) *sin (Om_y) ‘sin (A_y) +sin (Om_z) 
*sin (A_z) ‘sin (Om_x) *cos (A_x) *cos (Om_y) ) +sin (Psi) ‘sin (Phi) /cos (Psi) 

/ (cos (Phi) A 2+sin (Phi) A 2) *sin (Om_z) * (cos (Om_x) *cos (A_z) +sin (Om _x) 

*sin (A_x ) ‘sin (A_z) ) / (cos (Om_x) ‘cos (Ora_y) ‘cos (Om_z) +cos (Om_x) 

*sin (Om_y) *sin (A_y) *sin (Om_z) *cos (A_z) +sin (Om_y) *cos (A_y)'*sin (0 m_x) 
*sin (A_x) * cos (Om_z) -sin (Om_y) *cos (A_y) *sin (Om_x) *cos (A_x) *sin (Om_z) 
*cos (A_z) +sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin (Om_y) *sin (A_y) 
+sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A _x) *cos (Om_y) ) ) /S_y* (V_y-b_y) 

+ (sin (Om_x) * (sin (A_x) *sin (Om_y) *sin (A_y) +cos (A_x) *cos (Om_y) ) 

/ (cos (Om_x) *cos (Om_y) *cos (Om_z) +cos (Om_x) *sin (Om_y) *sin (A_y) 

*sin (Om_z) *cos (A_z) +sin (Om_y) *cos (A_y) *sin (Om_x) *sin (A_x) ‘cos (Om_z) 
-sin (Om_y) *cos (A_y) *sin (Om_x) *cos (A_x) *sin (Om_z) *cos (A_z) +sin (Om_z) 
*sin (A_z ) *sin (Om_x) *sin (A_x) *sin (Om_y) *sin (A_y) +sin (Om_z) *sin (A_z) 
*sin (Om_x) *cos (A_x) *cos (Om_y) ) +sin (Psi) *cos (Phi) /cos (Psi) / (cos (Phi) A 2 
+sin (Phi) A 2) *sin (Om_y) * (cos (Om_x) *sin (A_y) -sin (Om_x) *cos (A_x) 

*cos (A_y) ) / {cos (Om_x) *cos (Om_y) *cos (Om_z) +cos (Om_x) *sin (Om~y) 

*sin (A_y) *sin (Om_z) *cos (A_z) +sin (Om_y) *cos (A_y) *sin (Om_x) *sin (A_x) 
*cos (Om_z) -sin (Om_y) *cos (A_y) *sin (Om_x) *cos (A_x) *sin (Om_z) *cos (A_z) 
+sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin (Om_y) *sin (A_y) fsin (Om_z) 
*sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) +sin (Psi) *sin (Phi) /cos (Psi) 

/ (cos (Phi) A 2+sin (Phi) A 2) * (cos (Om_x) *cos (Om_y) +sin (Om_x) *sin (A_x) 

*sin (Om_y) *cos (A_y) ) / (cos (Om__x) *cos (Om_y) ‘cos (Om_z) +cos (Om_x) 

*sin (Om_y) *sin (A_y) *sin (Om_z) ‘cos (A_z) +sin (Om_y) *cos (A_y) *sin (Om_x) 
*sin (A_x) ‘cos (Om_z) -sin (Om_y) ‘cos (A_y) *sin (Om_x) ‘cos (A_x) *sin (Om_z) 
*cos (A_z) +sin (Ora_z) *sin (A_z) *sin(0m_x) *sin (A_x) *sin (Om_y) *sin (A_y) 
+sin (Om_z) *sin (A_z) *sin (Om_x) ‘cos (A_x) ‘cos (Om_y) ) ) /S_z* (V z-b z) ; 


f (2) = (cos (Phi) /cos (Psi) / (cos (Phi) A 2 + sin (Phi) A 2) *sin (Om_y) * (cos (A_y) 

* cos (Om_z) +sin (A_y) *sin (Om_z) *sin (A_z) ) / (cos (Om_x) ‘cos (Om_y) ‘cos (Om_z) 
+COS ( Om_x ) *sin (Om_y) *sin (A_y) ‘sin (Om_z) ‘cos (A_z) +sin (Om_y) ‘cos (A_y) 
‘sin (Om_x) ‘sin (A_x) ‘cos (Om_z) -sin (Om_y) ‘cos (A_y) ‘sin (Om_x) ‘cos (A_x) 
‘sin (Om_z) ‘cos (A__z) +sin (Om_z) ‘sin (A_z) ‘sin (Om_x) ‘sin (A _x) ‘sin (Om_y) 
‘sin (A_y) +sin (Om_z) ‘sin (A_z) ‘sin (Om_x) ‘cos (A_x) ‘cos (Om_y) ) +sin (Phi) 
/cos (Psi) / (cos (Phi ) A 2+sin (Phi) A 2) ‘sin (Om_z) * (-sin (Om_y) ‘cos (A_y) 

‘cos (A_z) +cos (Om_y) ‘sin (A_z) ) / (cos (Om_x) ‘cos (Om_y) ‘cos (Om z) 

+cos (Om_x) ‘sin (Om_y) ‘sin (A_y) ‘sin (Om_z) ‘cos (A_z) +sin (0m_y7*cos (A_y) 
‘sin (Om_x) ‘sin (A_x) ‘cos (Om_z) -sin (Om_y) ‘cos (A_y) ‘sin (Om_x) ‘cos (A_x) 
‘sin (Om_z) ‘cos (A_z) +sin (Om_z) ‘sin (A_z) ‘sin (Om_x) ‘sin (A_x) ‘sin (Om_y) 
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*sin (A_y) +sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) ) /S x 

* (V_x-b_x) + (cos (Phi) /cos (Psi) / (cos (Phi) A 2+sin ( Phi ) A 2 ) *7cos (OrrTx) 

*cos (Om_z) +sin (Om_x) *cos (A_x) *sin (Om_z) *sin (A_z) ) / (cos (Om x) 

* co s (Om_y) *cos (Om_z) +cos (Om_x) *sin (Om_y) *sin(A_y) *sin (Om_z) *cos (A_z) 
+sin (Om_y) *cos (A_y) *sin (Om_x) *sin (A_x) *cos (Om_z) -sin (Om_y) *cos (A_y) 
*sin(Om_x) *cos (A_x) *sin(Om_z) *cos (A_z) +sin (Om_z) *sin (A_z) *sin (Om_x) 
*sin (A _x) *sin (Om_y) *sin (A_y)+sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A_x) 
*cos ( Om_y ) ) -sin (Phi) /cos (Psi) /(cos (Phi) A 2+sin (Phi) A 2) *sin(Om_z)~ 

* (cos (Om_x) *cos (A_z) +sin (Om_x) *sin (A_x) *sin(A__z) ) / (cos (Om_x) 

*cos (Om_y) *cos (Om_z) +cos (Om_x) *sin (Om_y) *sin (A_y) *sin (Om_z) *cos (A_z) 
+sin (Om_y) *cos (A_y) *sin (Om_x) *sin (A_x) *cos (Om_z) -sin (Om_y) *cos (A_y) 
*sin (Om_x) *cos (A_x) *sin (Om_z) *cos (A_z) +sin (Om_z) *sin (A z) *sin (Om x) 
*sin (A_x) *sin (Om_y) *sin (A_y) +sin (Om_z) *sin (A_z) *sin (OnTx) *cos (A_x) 
*cos ( Om_y ) ) ) /S_y* (V_y-b_y) + ( -cos ( Phi ) /cos ( Psi ) / (cos (Phi) A 2 
+sin (Phi) A 2) *sin (Om_y) * (cos (Om_x) *sin (A_y) -sin (Om_x) *cos (A_x) 

*cos (A_y) ) / (cos (Om_x) *cos (Om_y) *cos (Om_z) +cos (Om_x) *sin (Om y) 

*sin (A_y) *sin (Om_z) *cos (A_z) +sin (Om_y) *cos (A_y) *sin (Om_x) *sin (A_x) 
*cos (Om_z) -sin (Om_y) *cos (A_y) *sin (Om_x) *cos (A_x) *sin (Om_z) *cos (A z) 
+sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin (Om_y) *sin (A_jy) +sin(OnTz) 
*sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) -sin (Phi) /cos (Psi) / (cos (Phi) A 2 
+sin (Phi) A 2) * (cos (Om_x) *cos (Om_y) +sin (Om_x) *sin (A_x) *sin (Om_y) 

*cos (A_y) ) / (cos (Om_x) *cos (Om__y) *cos (Om_z) +cos (Om_x) *sin (Om y) 

*sin (A_y) *sin (Om_z) *cos (A_z) +sin (Om_y) *cos (A_y) *sin (Om_x) *sin (A x) 
*cos (Om_z ) -sin (Om_y ) *cos (A_y ) *sin (Om_x ) *cos (A_x ) *sin (Om 2 ) *cos (A z ) 
+sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin (Om_y) *sin (A_y) +sin (Om^z) 
*sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) ) /S_z* (V z-b z); 


f ( 3) - ( sin ( Phi ) / (cos ( Phi ) A 2+sin ( Phi ) A 2 ) *sin (Om_y) * (cos (A_y) *cos (Om z ) 
+sin (A_y) *sin (Om_z) *sin (A_z) ) / (cos (Om_x) *cos (Om_y) *cos (Om_z) 

+cos ( Om_x ) *sin (Om_y) *sin (A_y) *sin (Om_z) *cos (A_z) +sin (Om_y7*cos (A_y) 
*sin (Om_x) *sin (A _x) *cos (Om_z) -sin (Om_y) *cos (A_y) *sin (Om x) *cos (A _ x) 
*sin (Om_z) *cos (A_z) +sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin(Om7y) 
*sin (A_y) +sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) -cos (Phi) 
/ (cos (Phi) A 2 + sin (Phi) A 2) *sin (Om_z) * (-sin (Om_y) *cos (A_y) *cos (A z) 

+ cos (Om_y) *sin (A_z) ) / (cos (Om_x) *cos (Om_y) *cos (Om_z) +cos (Om x) _ 

*sin (Om_y) *sin (A_y) *sin (Om_z) *cos (A_z) +sin (Om_y) *cos (A y) *sin (Om x) 
*sin (A_x) *cos (Om_z) -sin (Om_y) *cos (A_y) *sin (Om_x) *cos (A^x) *sin (Om _ z) 
*cos (A_z) +sin (Om_z) *sin(A_z) *sin (Om_x) *sin (A_x) *sin(Om_y) *sin(A_y) 
+sin (Om_z ) *sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) ) /S_x* (V x-b x7 
+ (sin (Phi) / (cos (Phi) A 2+sin (Phi) A 2) * (cos (Om_x) *cos (Om_z) +sln (Om x) 
*cos (A_x) *sin (Om_z) *sin(A_z) ) / (cos (Om_x) *cos (Om_y) *cos (Om_z) 

+cos (Om_x) *sin (Om_y) *sin (A_y) *sin (Om_z)*cos (A_z) +sin (Om_y7*cos (A y) 
*sin (Om_x) *sin (A_x) *cos (Om_z) -sin (Om_y) *cos (A_y) *sin (0m7x) *cos (A _ x) 
*sin (Om_z) *cos (A_z) +sin (Om_z) *sin (A_z) *sin (0m_x) *sin (A_x) *sin (Om _ y) 
*sin (A_y ) +sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A _x) *cos (Om_y) ) +cos (Phi) 
/ (cos (Phi) A 2 + sin (Phi) A 2) *sin (Om_z) * (cos (Om_x) *cos (A_z) +sin (Om x) 
*sin (A_x) *sin (A_z) ) / (cos (Om_x) *cos (Om_y) *cos (Om_z) +cos (Om x) 

*sin (Om_y) *sin (A_y) *sin (Om_z) *cos (A_z) +sin (Om_y) *cos (A y) *sin (Om x) 
*sin (A_x) * cos (Om_z) -sin (Om_y) *cos (A_y) *sin (Om_x) *cos {A~x) *sin (Om _ z) 
*cos (A_z) +sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin (Om_y) *sin(A y) 
+sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) ) /S_y* (V y-b y7 
+ (-sin (Phi) / (cos (Phi) A 2+sin (Phi) A 2) *sin (Om_y) * (cos (Om_x) *sin (A y) 

-sin (Om_x) *cos (A_x) *cos (A_y) ) / (cos (Om_x) *cos (Om_y) *cos (Om z) _ 

+cos ( Om_x ) *sin (Om_y) *sin (A_y) *sin (Om_z) *cos (A_z) +sin (0m_y7*cos (Ay) 
*sin (Om_x) *sin (A_x) *cos (Om_z ) -sin (Om_y) *cos (A_y) *sin (Om x) *cos (A~x) 
*sin (Om_z) *cos (A_z) +sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin (Om^y) 
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*sin (A_y) +sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A__x) *cos (Om_y) ) 

+cos (Phi) / (cos (Phi) A 2+sin (Phi) A 2) * (cos (Om_x) *cos (Om_y) +sin(Om_x) 
*sin (A_x) *sin (Om_y) *cos (A_y) ) / (cos (Om_x) *cos (Om_y) *cos (Om_z) 

+cos (Om_x) *sin (Om_y) *sin (A_y) *sin (Om_z) *cos (A_z) +sin (Om_y) *cos (A_y) 
*sin (Om_x) *sin (A_x) *cos (Om_z) -sin (Om_y) *cos (A_y) *sin (Om_x) *cos (A_x) 
*sin (Om_z) *cos (A_z) +sin (Om_z) *sin (A_z) *sin (Om_x) *sin (A_x) *sin (Om_y) 
*sin (A_y) +sin (Om_z) *sin (A_z) *sin (Om_x) *cos (A_x) *cos (Om_y) ) ) /S_z 
* (V z-b z ) ; 


return; 
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